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Most  chemical  reactions  take  place  in  solution  and,  in  order  to  understand 
the  molecular  processes  that  take  place  during  reaction,  it  is  necessary  to  include 
the  effect  of  the  surrounding  medium  on  the  reacting  species.  As  the  solvent 
is  usually  in  molar  quantities,  it  is  impossible  to  study  an  Avagadro’s  number  of 
molecules  using  the  principles  of  quantum  mechanics.  To  overcome  this  difficulty, 
reasonable  models  that  adequately  describe  the  solvent  have  to  be  constructed. 

Three  such  models  are  presented  in  this  thesis.  The  first  describes  the  solvent 
as  a continuum  with  the  solute  molecule  centered  in  a spherical  cavity  embedded 
in  the  solvent.  This  model  does  not  assign  a structure  to  the  molecules  in  the 
solvent.  The  second  model  builds  a cluster  of  solvent  molecules  around  the  solute 
using  the  method  of  simulated  annealing,  a method  based  on  the  Metropolis  Monte 
Carlo  algortihm.  The  specific  nature  of  the  interactions  between  the  solute  and 
solvent  molecules  is  taken  into  consideration  in  this  method. 


VI 


The  third  model  uses  the  second  model  to  build  a cluster  of  solvent  molecules 
around  the  solute  and  this  entire  solute-solvent  cluster  is  then  enclosed  in  a 
spherical  cavity  and  embedded  in  the  solvent. 

Solvent  induced  spectral  shifts  in  solute  molecules  were  evaluated  using 
the  first  and  third  models  and  the  results  obtained  compared  with  experimental 
shifts.  The  method  used  to  calculate  the  shifts  was  the  Intermediate  Neglect 
of  Differential  Overlap  which  has  been  parameterized  for  the  study  of  spectral 
transitions. 

Not  all  molecules  can  be  described,  even  approximately,  by  a spherical  cavity 
model.  An  algorithm  is  presented  that  defines  a more  realistic  shape  for  this  cavity 
as  a set  of  interlocking  spheres  centered  on  the  atoms  of  the  solute.  The  results 
of  the  application  of  this  method  to  such  processes  as  molecular  dissociation  are 
also  presented. 


CHAPTER  1 
INTRODUCTION 


Solvation  is  the  process  in  which  one  system,  usually  at  a lower  concentration, 
interacts  with  another  system,  which  is  in  larger  amounts.  The  first  system  is 
referred  to  as  the  solute  and  the  latter  is  termed  the  solvent.  The  interaction 
between  the  two  systems  can  affect  the  structure  of  the  constituent  molecules. 
Consequently,  the  chemical  and  physical  properties  of  these  molecules  show  a 
change  from  their  corresponding  values  in  the  gas  phase  or  vacuum. 

An  exact  numerical  study  of  solvent  effects,  based  on  the  principles  of 
quantum  mechanics  is  impossible,  even  after  the  advent  of  supercomputers  as 
the  number  of  electrons  and  nuclei  which  constitute  the  interacting  molecules 
are  of  the  order  of  Avogadro’s  number.  Even  an  approximate  theory  for  such 
a number  of  particles  is  questionable,  as  the  simplifications  would  have  to  be 
drastic,  thereby  making  the  interpretation  of  the  results  difficult 

In  order  to  overcome  this  difficulty,  models  that  adequately  represent  the 
solvent  have  to  be  constructed  and  then  incorporated  into  the  quantum  mechanical 
description  of  the  solute  system  (1).  The  validity  of  these  models  is  usually 
checked  by  evaluating  thermodynamic  properties  like  the  specific  heats  and 
isothermal  compressibility,  and  comparing  these  results  with  experiment. 
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Solvent  models  can  be  broadly  classified  into  one  of  the  following  categories: 

1.  Continuum  models 

2.  Discrete  models. 

A synthesis  of  the  two,  referred  to  as  the  discrete -continuum  approach, 
incorporates  the  positive  aspects  of  each  model  for  a better  description  of  the 
solvent. 

In  the  continuum  model,  the  solute  molecule  is  represented  by  a point  dipole 
centered  in  a cavity  which  is  embedded  in  the  solvent.  The  shape  of  this  cavity  can 
be  chosen  beforehand  and  usually  is  taken  to  be  a sphere,  ellipsoid  or  a rod-like 
tube.  The  solvent  is  described  as  a continuous  medium  with  a specific  dielectric 
constant  (e)  and  an  associated  refractive  index  (n).  Clearly,  this  description  does 
not  take  the  structure  of  the  solvent  molecules  into  account. 

The  specific  nature  of  the  solute- solvent  interactions  is  important,  particularly 
in  the  case  of  polar  species.  The  dipolar  interactions  are  most  important  in  ordering 
the  solvent  molecules  around  the  solute  to  form  solvation  shells.  The  properties 
of  the  molecules  in  these  shells  are  quite  different  from  those  in  the  bulk  areas 
of  the  solvent.  By  using  the  same  value  for  the  dielectric  constant,  the  specific 
nature  of  the  solute  solvent  interactions  is  neglected. 

Beveridge  and  Schnuelle  (2)  attempted  to  include  a solvation  shell  within  the 
continuum  model  by  centering  the  point  dipole  inside  two  concentric  spheres.  The 
value  of  the  dielectric  constant  for  the  solvent  between  the  two  spheres  is  differ- 
ent from  that  of  the  bulk  solvent  which  is  outside  the  second  sphere.  The  major 
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shortcomings  in  this  model  are  the  choice  of  values  for  the  dielectric  constants 
and  the  radii  defining  the  spheres,  which  are  arbitrary. 

Surrounding  the  solute  with  a suitable  number  of  solvent  molecules  is  the  basis 
of  the  discrete  model.  The  number  of  such  molecules  is  chosen  based  on  intuition 
and  prior  experience.  Studies  on  such  solute-solvent  complexes  explicitly  include 
the  anisotropic  nature  of  the  associated  interactions.  However,  the  disadvantage 
in  using  this  scheme  is  the  choice  of  an  adequate  number  of  solvent  molecules 
that  accurately  describe  the  solvent.  For  typical  dipolar  cases,  this  number  can  be 
in  the  thousands  which  makes  a meaningful  numerical  treatment  impossible. 

The  solvent  model  of  Noell  and  Morokuma  (3)  is  obtained  from  good  ab-initio 
calculations  on  supermolecular  complexes.  The  atoms  in  the  solvent  molecules 
are  then  replaced  by  fractional  point  charges  which  are  chosen  to  reproduce  the 
interaction  energy  calculated  for  the  complex.  The  drawbacks  in  this  approach 
are  that  a prior  knowledge  of  the  solvent  structure  is  required  and  the  effects  of 
the  bulk  solvent  are  now  neglected. 

A natural  extension  of  the  two  models  described  above  is  the  discrete- 
continuum  approach.  Based  on  a suitable  theoretical  formulation,  the  solute  is 
surrounded  by  a predetermined  number  of  solvent  molecules  and  the  resulting 
system,  which  now  consists  of  explicit  solvation  shells,  is  then  embedded  in  the 
fluid  continuum.  The  quantum  mechanical  study  of  the  solvated  complex  within 
the  discrete  model  can  be  achieved  using  the  following  methods: 
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1.  Variation  scheme. 

2.  Perturbational  treatment. 

In  the  former  case,  the  solvated  complex  is  treated  as  one  molecule  referred 
to  as  a supermolecule.  The  interaction  energy  AE  of  the  solute  molecule  A with 
the  N surrounding  solvent  molecules  Bi,  {i=l,2,  ...,  N}  is  calculated  as 


The  first  term  on  the  right  hand  side  of  the  above  equation  is  the  energy 
of  the  supermolecule.  Ea  and  Eei  are  the  energies  of  the  solute  and  the  i* 
solvent  molecule  respectively.  The  last  term  in  the  equation  is  the  solvent  — 
solvent  interaction  energy.  These  energies  are  usually  calculated  in  the  Hartree- 
Fock  scheme  using  a finite  basis  set  Clearly,  the  supermolecule  approach  is 
practical  only  when  the  size  of  the  solute  and  the  number  of  solvent  molecules 
are  very  small,  and,  even  in  such  cases,  the  choice  of  the  basis  set  is  very 
important  in  obtaining  meaningful  results  (4).  The  numerical  accuracy  required 
in  determining  each  of  the  terms  on  the  right  hand  side  of  the  equation  (1.1)  has 
to  be  extraordinarily  high  as  the  interaction  energy  is  obtained  by  subtraction. 

If  the  solvent  induced  structural  changes  in  the  solute  are  considered  to  be 
small,  a perturbational  scheme  seems  to  be  more  suitable  in  studying  intermolec- 
ular  interactions.  Using  the  theory  of  weak  interactions,  the  energy  of  interaction 
AE  can  be  expanded  in  a multipolar  series  expressed  as  a function  of  the  inverse 
of  the  distance  separating  the  systems  (5).  An  advantage  in  such  a formalism 
is  that  each  term  in  the  series  has  a physical  meaning,  such  as  the  coulombic, 
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induced  energy  contributions,  etc.  This  generalized  multipolar  expansion  is  an  in- 
finite series  and  for  reasons  of  computational  efficiency  has  to  be  truncated  after 
a suitable  number  of  terms.  The  coefficients  of  each  of  the  remaining  terms  are 
determined  by  fitting  the  expansion  to  experimental  thermodynamical  quantities 
(6)  or  values  of  the  interaction  energy  obtained  for  model  systems  from  good 
ab-initio  calculations  (7). 

After  the  multipolar  expansion  for  AE  has  been  defined,  the  interaction  energy 
for  the  solute-solvent  cluster  can  be  evaluated.  The  number  of  solvent  molecules 
in  such  a cluster  can  be  determined  by  calculating  the  interaction  energy  per 
solvent  molecule,  AEi  for  an  increasing  number  of  solvent  molecules  (m)  around 
the  solute.  The  value  of  m at  which  AEj  converges,  within  a predetermined 
value,  is  chosen  to  be  the  number  of  solvent  molecules  in  the  solvation  shell.  The 
solvated  complex  is  then  centered  in  a cavity  of  suitable  dimensions  and  emersed 
in  the  continuum  fluid.  In  this  way,  the  anisotropy  in  the  solvent  shells  and  the 
interactions  due  to  the  bulk  fluid  can  be  accounted  for. 

Shifts  in  the  electronic  spectra  of  molecules  due  to  solvation  is  the  main  area 
of  research  presented  in  this  dissertation.  Certain  dyes  can  have  different  colors 
depending  on  the  nature  of  the  solvent  they  are  in.  An  example  is  4-nitro  - 4'  - 
(N,N  diethyl)  diazabenzene  (Figure  1.1)  which  is  yellow-orange  in  cyclohexane 
and  deep  red  in  ethanol. 

Empirical  rules  have  been  proposed  (8,9)  to  calculate  these  shifts.  Although 
these  rules  are  helpful  in  predicting  the  positions  of  spectral  peaks  in  various  sol- 
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vents  to  within  5 to  10  nm  of  experiment,  they  offer  no  insight  into  the  molecular 
processes  that  cause  these  shifts.  An  understanding  of  these  interactions  should 
be  based  on  good  microscopic  concepts. 


Figure  1.1  Structure  of  4-nitro  - 4'-  (N,N  diethyl)  diazabenzene. 

In  order  to  study  the  origin  of  such  spectral  shifts,  the  three  solvent  models 
described  above  were  used.  In  each  case  the  effective  nature  of  the  surroundings 
was  incorporated  into  the  Hamiltonian  of  the  solute.  The  solute  Hamiltonian  is 
obtained  from  the  Intermediate  Neglect  of  Differential  Overlap  (INDO)  approx- 
imation (10)  to  the  Hartree-Fock  self-consistent  field  equations  (11)  as  parame- 
terized by  Ridley  and  Zemer  (12)  and  Bacon  and  Zemer  (13),  but  the  methods 
developed  are  not  limited  to  this  model. 
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Chapter  two  of  this  thesis  introduces  the  partitioning  of  the  solute-solvent 
system.  The  solvent  is  treated  as  a continuum  and  the  mathematical  formulation 
for  the  calculations  derived.  The  values  calculated  for  the  shifts  in  the  spectral 
peaks  for  a series  of  molecules  in  various  solvents  are  tabulated. 

Obtaining  the  discrete  solvent  strucure  is  the  main  focus  of  the  theory  outlined 
in  Chapter  3.  The  theory  of  the  Simulated  Annealing  method  (14)  and  the  results 
obtained  for  the  solvation  of  neutral  molecules  and  ions  are  presented. 

The  structures  obtained  for  the  solvent  clusters  using  the  method  described  in 
chapter  three  are  replaced  by  point  charges  having  the  same  values  as  those  used 
in  the  simulation.  The  spectral  shifts  in  such  clusters  are  calculated  and  compared 
with  the  corresponding  experimental  values  obtained  from  laser  jet  studies.  These 
solute-solvent  clusters  are  then  embedded  in  the  continuum  solvent.  The  theory 
used  to  calculate  the  shifts  within  the  discrete-continuum  model  along  with  the 
results  obtained  are  presented  in  Chapter  4. 

An  algortihm  which  defines  a more  realistic  shape  of  the  molecular  cavity  in 
the  solvent  as  a set  of  interlocking  van  der  Waals  spheres  centered  on  the  atoms  of 
the  solute  is  presented  in  Chapter  5.  The  effect  of  the  image  charges  induced  by 
the  solute  on  the  surface  of  the  cavity  is  incorporated  into  the  solute  Hamiltonian. 
Processes  like  molecular  dissociation  are  studied  and  the  results  tabulated. 

The  last  chapter  gives  the  conclusions  reached  on  the  use  of  the  methods 
presented  in  this  thesis.  The  advantages  and  shortcomings  in  the  methods  and  the 
parameterization  are  outlined,  along  with  possible  extensions  for  the  future. 


CHAPTER  2 

SPECTRAL  SHIFTS  USING  THE  CONTINUUM 
MODEL  FOR  THE  SOLVENT 

2.1  Introduction 

The  configuration  space  spanned  by  the  solute  and  the  solvent  can  be  parti- 
tioned into  three  subspaces  — one  that  is  described  by  the  electrons  and  nucleii 
of  the  solute  S,  the  second  subspace  D,  that  is  spanned  by  the  discrete  solvent 
molecules  immediately  surrounding  the  solute,  and  the  continuum  solvent  space 
C,  described  by  the  dielectric  constant  (e)  and  the  refractive  index  (n).  As  the  sol- 
vent induced  changes  in  the  properties  of  the  solute  are  considered  to  be  small,  the 
effect  of  the  solvent  can  be  written  as  a perturbational  expansion  in  the  parameter 
A.  The  Hamiltonian  Hs  for  the  solvated  complex  can  now  be  written  as 

H'S  = HS  + X H(J1d  + X + X 3H^1d  + ... 

+AJ4-C  + x2hS-C  + xlHS-C  + ...  + A 2.1 

where  Hu,  I,J  e { S,  D,  C}  is  the  Hamiltonian  for  the  interaction  between  the 
subspaces  I and  J and  Hs  describes  the  quantum  mechanics  of  the  solute. 

A represents  the  interaction  between  the  continuum  and  discrete  subspaces,  and 
each  of  these  with  itself. 

In  this  chapter,  the  solvent  induced  shifts  in  the  electronic  spectra  of  a series 
of  molecules  are  calculated.  The  solvent  is  treated  as  a continuum;  that  is,  the 
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coupling  between  the  solute  space  and  discrete  solvent  space  is  neglected.  As 
the  time  elapsed  during  electronic  transitions  is  of  the  order  of  10  ~18  seconds, 
changes  in  the  solvent  structure  that  effect  the  interactions  within  the  solvent  can 
be  neglected.  Thus,  A is  assumed  to  be  the  same  for  the  ground  and  excited 
states  of  the  solute  and  need  not  be  considered  explicitly  in  the  formulation.  The 
approximate  Hamiltonian  H”s  for  the  solvated  complex  is  now  given  by 

ffS  = HS  + A H^c  + A *n{P_c  + A *Hf_c  + ...  2.2 

Expressions  for  the  solute-solvent  interaction  are  derived  based  on  the  theory  of 
weak  interactions  and  the  energy  of  interaction  is  then  parameterized  using  the 
reaction  field  theory  of  Onsager  (15).  These  derivations  are  given  in  sections  2.2 
and  2.3  respectively. 

The  Intermediate  Neglect  of  Differential  Overlap  (INDO)  method  (10)  has 
been  particularly  useful  in  the  prediction  of  electronic  spectra.  This  method,  as 
parameterised  by  Ridley  and  Zemer  (12)  and  Bacon  and  Zemer  (13)  is  used  to 
describe  the  solute  molecule.  The  formulation  of  the  method  is  well  documented 
and  will  not  be  repeated  here. 

The  incorporation  of  the  effective  solvent  Hamiltonian  and  derivation  of 
the  wavefunction  of  the  solute  (in  the  solvent)  is  presented  in  section  2.4, 
the  restrictions  on  the  wavefunction  during  the  variation  are  the  orthonormality 
requirement  and  the  reaction  field  energy  derived  in  section  2.3. 
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The  Onsager  model  treats  the  solute  molecule  as  a point  dipole  centered  in  a 
cavity.  The  shape  of  the  cavity  chosen  here  is  a sphere  and  the  methods  used  to 
evaluate  the  radius  of  this  cavity  are  discussed  in  section  2.5 

Section  2.6  gives  the  calculated  and  experimental  spectral  shifts  for  a series 
of  molecules  in  various  solvents.  The  final  section  gives  the  conclusions  reached 
from  this  study  and  criticisims  of  the  continuum  model. 

2.2  Derivation  of  Equation  for  the  Interaction  Energy  AE 


When  the  distance  separating  two  interacting  systems  A and  B is  larger  than 
their  individual  dimensions,  the  supermolecular  system  can  be  studied  using  the 
theory  of  weak  interactions.  Neglecting  the  permutation  of  electrons  between 
the  systems  is  justified  on  the  basis  of  the  large  intersystem  distance  and  the 
wavefunction  of  the  composite  system  can  be  approximated  as  the  simple 
product  of  the  wavefunctions  of  the  individual  systems  'I' a and  as 

* = {Aa^a){Ab^b) 

where  Aa  and  Ab  are  the  antisymmetrizers  for  the  electrons  in  systems  A and 
B respectively. 

In  the  following  derivation,  the  solute  whose  spectrum  is  to  be  studied  is 
denoted  by  A,  while  the  N solvent  molecules  surrounding  it,  which  comprise 
system  B,  will  be  referred  to  as  Bi,  B2,  . . .,  Bn.  Let  a and  0 refer  to  the  nuclei 
and  a,  b to  the  electrons  of  systems  A and  B^,  { k = 1,  2,  . . .,  N } respectively. 
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Further,  let  these  particles  have  coordinates  ra,  r/9,  ra  and  rb  relative  to  origins 
Oa  and  C>Bk  with  the  distance  0AC>Bk  being  equal  to  Rk- 

The  Hamiltonian  H for  the  supermolecular  system  ABk  can  be  written  as 


H = Ha  + HBk  + HABk 


2.3 


where  Hi,  ( I = A,  Bk  ) and  HABk  are  the  Hamiltonians  for  system  I and  the 
interaction  between  A and  Bk  respectively.  HABk  can  be  written  in  the  usual 
manner  as 


HABk 


ST'  ZgZp  Zj  Zi  y>  1 24 

ra/3  ri'  i«inri'  rab 

a<p  l=oi  t=o  l=p  t=a  a<b 


Z\  (I  = a,  /?)  is  the  charge  on  the  nucleus  I.  The  distance  rMN  between  particles 
M and  N,  where  M,N  can  be  nucleii  or  electrons  is  given  by 

nn  , n(zN  ~ ZM) 

rMN  — + 2 — 


2 2 2 
(xn-xm)  +(yN-VM ) +(zn~zm)  -.3-  ?4- 

If  the  two  systems  A and  Bk  are  electronically  neutral,  equation  (2.5)  substituted 
in  equation  (2.4)  yields,  after  manipulation,  a leading  term  in  R-3  as  the  R-1 
and  R-2  terms  disappear.  Defining  the  dipole  moment  operators  MA  and  MBk  of 
systems  A and  Bk  as 
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MA  = ^2  Z<*r<*  ~ y,  ra 


a£A  a£A 


MBk  = ^2  zPrp  - 5Z rb 


0eBk  beBk 


HABk  can  be  written  as 


HABk  = R~3MAeABkMBk  + O (JT4) 


2.6 


where  0ABk  is  the  orientation  tensor  between  A and  Bk.  The  Hamiltonian  Hab  for 
the  solute  molecule  in  the  presence  of  the  N solvent  molecules  is  now  given  by 


Equations  (2.6)  and  (2.7)  indicate  that  the  leading  term  of  the  interaction  is  the 
dipole -dipole  term. 

Evaluation  of  the  interaction  energy  AE  by  treating  the  interaction  Hamil- 
tonian Hab  as  a perturbation  requires  the  definition  of  a reference  wavefunction 
for  the  solute- solvent  complex.  To  evaluate  the  second  order  interaction  en- 
ergy, excited  state  wavefunctions  for  the  solvated  complex  need  to  be  defined. 
From  equation  (2.7),  it  is  clear  that  Hab  involves  a summation  over  the  solvent 
molecules  Bk,  one  molecule  at  a time.  It  is  therefore  necessary  to  consider  only 
those  excited  state  wavefunctions  in  which  the  solute  molecule  and/or  one  solvent 
molecule  is  excited. 

Assuming  that  the  solute  and  solvent  molecules  are  far  apart,  the  wavefunc- 


N 


2.7 


tions  for  the  solvated  solute  system  in  the  weak  interaction  limit  are 
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*0  — V’^l/,5lV,52---V,5Jt-"'05Ar  2.8a 

^1  = {V’yll/,51,/,B2-"V,5it---V’B7v}  2.8b 

^2  = {V,^V’5l^,fi2-"V,Bit---'05Af}  2.8c 

^3  = 2.8d 


In  equations  (2.8  a-d),  ?/>%  and  V,0Bk  are  the  wavefunctions  for  the  solute  molecule 
and  the  k*  solvent  molecule  in  the  ground  state  (denoted  0).  The  excited  state 
wavefunctions  for  these  molecules  are  denoted  by  t/>ca  and  V>bBk  respectively. 

= { },  1=1,  2,  3,  represents  the  set  of  all  possible  inferred  excitations. 

The  first  order  change  in  the  energy  E(1)  due  to  the  solute-solvent  interaction 
Hab  is  given  by 

£(1)  =<  <b0\HAB\<Ho  > 


Substituting  equation  (2.8a)  in  the  above  equation  and  using  the  definitions  of 
the  dipole  moment  operators  MA  and  MBk  given  before 

P=1 

where  /z-4  (I=A,  B ; j = 0,  e,  b ) is  the  dipole  moment  of  system  I in  state  j. 
E(2),  the  second  order  change  in  the  energy,  is  given  by 


2.9 


£(2)  = - E 


/=! 


< <Hq\Hab\^,  > 
(E,  - Eq) 


2.10 
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where  the  wavefunctions  $/,  / = 1,2,  3,  have  been  defined  in  equations  (2.8b  — 
d).  Defining  the  polarisability  tensor  o°a  of  A in  state  0 as 

< > |2 


aA  = 


e/0 


(E\  ~ E°a ) 


2.11 


the  second  order  contribution  to  the  energy  from  wavefunction  'I'  \ is  given  by 

i N 

171(2)  ^ \ ' p— 3 p— 3 0 f\ABp  0 r\ABq  0 7 17 

E 1 L RV  Ro  Pq b ,pu  qfiB,g  ZAl 


r,9=i 


Similarly,  the  contribution  E2(2)  from  $2  is 


N 


E22)  = -5  E 


2.13 


P=1 


where  the  polarisability  tensor  of  B is  defined  analogous  to  equation  (2.11). 
E3 (2)  which  is  defined  as 


p(2)  _ 1 < ^q\HAB\^3  > 

3 (£3  - E0 ) 


is  the  energy  due  to  the  interaction  of  the  excited  states  of  the  solute  and  solvent 
molecules  and  this  term  is  the  dispersion  energy  contribution  to  the  second  order 
energy.  In  the  method  presented  here,  the  explicit  form  of  this  term  will  be 
neglected  and  will  henceforth  be  denoted  as  Dg. 

The  change  in  energy  AE  of  the  solute,  to  second  order,  can  now  be  written  as 


A E = + E j2)  + E{2)  + Dg 


2.14 


where  the  terms  on  the  right  hand  side  of  the  above  equation  have  been  defined 
in  equations  (2.9),  (2.12)  and  (2.13). 
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Each  of  the  terms  in  equation  (2.14)  is  calculated  as  a sum  over  the  N solvent 
molecules.  When  the  number  of  such  molecules  is  of  the  order  of  Avagadro’s 
number,  such  summations  are  not  practical  and  have  to  be  replaced  by  statistical 
mechanical  averages  over  all  the  configurations  of  the  solvent  molecules,  leading 
then  to  macroscopic  parameters. 

2.3  Parameterization  of  the  Expression  for  AE 


The  reaction  field  theory  of  Onsager  (15)  treats  the  solvent  as  a continuum 
and  the  inherent  assumption  is  that  the  solvent  is  isotropic.  In  other  words,  no 
molecular  structure  is  assigned  to  the  solvent.  Thus,  the  polarisability  tensors  in 
the  above  equations  can  be  replaced  by  their  corresponding  scalar  values.  As 
the  solvent  lacks  a structure,  the  summations  in  these  equations,  which  are  over 
the  number  of  solvent  molecules,  are  replaced  by  statistical  averages  in  terms  of 
macroscopic  properties  like  the  dielectric  constant  (e)  and  the  refractive  index  (n). 

The  solute  is  treated  as  a polarizable  dipole  centered  in  a cavity  which  is 
embedded  in  the  continuum.  For  purposes  of  simplicity,  the  cavity  is  assumed  to 
be  a sphere  of  radius  'a'.  The  solute  dipole  fi° a creates  an  image  dipole  in  the 
opposite  direction  in  the  solvent.  This  image  dipole,  in  turn,  augments  the  solute 
dipole  and  creates  a uniform  field  called  the  reaction  field,  R,  which  enhances  the 
already  existing  field  in  the  cavity. 

The  reaction  field  R is  given  by  (15) 

2 (D'  — 1)  o 

a3(2D'  + l)^ 


2.15 
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where  D'  is  that  part  of  the  static  dielectric  constant  which  takes  the  interactions 
between  the  solute  dipole  and  the  permanent  dipoles  of  the  solvent  molecules  into 
account.  The  derivation  of  this  equation  is  given  in  the  Appendix. 

Let 

r 2 (D1  — 1) 

1 a3  (2D1  + 1) 

For  a polarizable  solute  dipole  a with  isotropic  polarizability  a,  the  reaction 
field  R is  given  by 

R = f(n°A  + aR)  2.16 


Approximating  a by  a3/2  (16),  the  reaction  field  for  a polarizable  solute  is 
now  given  by 


D'  - 1 
D'  + 2 


A 


2.17 


The  static  dielectric  constant  e of  the  solvent  is  comprised  of  two  parts,  one 
due  to  the  permanent  dipole  orientation  D’  and  the  other  D",  due  to  the  electronic 
polarizability  of  the  solvent  (17,18).  The  latter  contribution  can  be  evaluated  using 
the  Maxwell  relation  D"  = n2,  where  n is  the  refractive  index  of  the  solvent. 

The  dielectric  term  in  the  above  equation  can  be  calculated  using  the  following 
relation: 

D'  - 1 e - 1 n2  - 1 
D'  + 2_e  + 2-n2  + 2 


2.18 
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Substituting  equation  (2.18)  in  (2.17)  yields 

„ 2 ft- 1 r?  - 1 ] 0 

R~  a3\t  + 2~  n2  + 2rA 


2.19 


Analysis  of  equation  (2.9)  shows  that  the  first  order  term  E(1)  is  due  to  the 
coupling  of  the  solute  dipole  with  the  permanent  dipoles  of  the  solvent  which 
is  exactly  the  reaction  field  R.  Therefore,  taking  the  statistical  average  over  the 
solvent  molecules,  the  term  E(1)  in  equation  (2.14)  can  be  written  as 

£{1)  = -h°aR  2.20 

where  R is  given  by  equation  (2.19). 

In  a similar  manner,  the  second  term  on  the  right  hand  side  of  equation 
(2.14)  is  the  coupling  of  the  field  due  to  the  permanent  solvent  dipoles  with  the 
polarizability  tensor  of  the  solute.  This  term  can  be  written  as  (19) 

e\2)  = -5  E WA^Vb,,,  2-21 

P>9=1 

The  field  at  the  p*  solvent  molecule  due  to  the  permanent  dipole  of  the  solute  is 
represented  by  the  third  term  of  equation  (2.14).  As  in  the  previous  case, 

42)  = -lT.Rp‘a°B,p(A)2  2.22 

p=l 

where  a°p  is  the  static  polarizability  of  the  pth  solvent  molecule.  Equation 
(2.22)  represents  the  induction  energy  due  to  the  polarization  of  the  solvent  by 
the  permanent  solute  dipole. 
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Averaging  the  terms  in  equation  (2.22)  over  the  molecular  orientations,  the 
reaction  field  r at  the  pth  solvent  molecule  due  to  the  induction  effects  of  the 
permanent  dipole  of  the  solute  molecule  is  given  by 


r = 


2.23 


where  n is  the  refractive  index  of  the  solvent. 

Thus,  AE  is  now  given  by 

AE  = -&R  - \ Y,  R?R-q3A,AA,r 

p>?=  1 


2.24 


The  time  scales  for  electronic  excitations  are  of  the  order  of  10  -18  seconds, 
about  a million  times  faster  than  the  vibrational  motion  of  the  nuclei.  When  the 
solute  molecule  is  excited,  its  dipole  is  altered  and  the  dipoles  of  the  solvent 
molecules  will  change  in  order  to  accommodate  this  effect.  The  reorientation  of 
the  permanent  part  of  the  solvent  dipole  can  be  neglected  as  the  nuclei  are  assumed 
to  be  fixed  during  the  electronic  transition.  Thus,  the  excitations  considered  obey 
the  Frank-Condon  principle.  However,  the  electronic  part  of  the  dipole  moment 
will  change,  and  this  effect  is  taken  into  consideration  by  using  the  reaction  field 
term  which  depends  on  the  refractive  index  of  the  solvent,  namely  equation  (2.23). 
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For  the  excited  state  of  the  solute  with  a dipole  moment  //a,  the  change  in 
energy  AEex  due  to  solvation  is  given  by 


where  De  is  the  dispersion  energy  term  for  the  excited  state  of  the  solute.  The  shift 
in  the  spectral  peak  Au  due  to  solvation  can  now  be  calculated  from  equations 
(2.24)  and  (2.25)  as 

Au  = A E — A Eex  2.26 


2.4  Variation  of  a Trial  Solute  Wavefunction  'l/'  in  the  Presence  of  the 

Continuum  Solvent 


Let 


and 


2 f 2(e  — 1)  n2-l] 

5 “ a3  \ e + 2 _ n2  -f  2 J 


2.27 


H*  = HS-  i/i T.g.fi  2.28 

where  H*  and  Hs  are  the  Hamiltonian  operators  for  the  solute  molecule  in  the 
solvent  and  gas  phase  respectively  and  /x  is  the  total  dipole  moment  operator. 
The  superscript  T denotes  the  transpose. 
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Define  a function  L according  to 

L = E-  A [<  > -1] 


+ |{<  > -9.  < tf'lAl*'  > ~9  {A)2} 

where  A and  7 are  constants  to  be  determined  and  is  the  normalized  trial 
wavefunction  for  the  solute  in  the  solvent.  /z° a is  the  expectation  value  of  the 
operator  /i  for  the  wavefunction  E is  the  energy  of  the  solute  molecule  in  the 
presence  of  the  continuum  solvent  and  is  given  by 

E=<  tf'lfTI#'  > 


Using  variational  theory 

6L  =<  $¥'|1T|¥'  > + < > 


-A  [<  6V’\*'  > + < *'{6*'  >] 


+^{<  > .g.  < > 2.29 


+ < > .g.  < > 


+ < > .g.  < > 
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+ < ^'\flT\S^'  > .g.  < *'1/21*'  >} 

At  a stationary  point  SL  = 0.  As  /z  is  a hermitian  operator, 

< >=<  tf'|/z|1'/  > 

Collecting  terms,  equation  (2.29)  reduces  to 

0 =<  — X + ig  M)2|^'> 


+ <V'\H*  - \ + jg(n°A)2\W  > 2.30 

As  the  variation  in  is  complex  and  arbitrary, 

<6*'\H*-\  + ig  (^)2|^'  >=  0 

and 

<¥'|ir-A  + 7<7  (^)2|^'>=0 


Identifying  A with  Es,  the  energy  of  the  solute  in  the  gas  phase,  and  7 = 1/2, 
*n*'>=  2.31 

Thus,  if  the  Hamiltonian  of  the  solute  is  defined  as  in  equation  (2.28),  the  energy 

E of  the  solute  in  the  presence  of  the  continuum  solvent  is  given  by 

1 2 


E = E„ 


MV 


e — 1 

7+2 


n~  — 1 
n2  + 2 


2.32 
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Clearly,  the  solvent  induced  energy  change  in  the  solute  is  due  to  the  reaction 
fields,  one  arising  from  the  interaction  of  the  permanent  dipoles  and  the  other 
due  to  the  coupling  of  the  solute  dipole  with  the  polarizability  of  the  solvent. 
Therefore,  the  corrections  to  the  gas  phase  energy  of  the  solute  include  corrections 
through  first  order  and  part  of  the  second  order  terms. 

Defining  the  Hamiltonian  as  in  equation  (2.28),  the  method  presented  above 
was  incorporated  into  the  INDO  formulation  of  Ridley  and  Zemer  (12)  and  Bacon 
and  Zemer  (13).  The  ground  state  wavefunction  I <f>'  > for  the  solute,  having 
2N  valence  electrons,  in  the  solvent  was  described  by  a single  determinant  in  the 
molecular  spin  orbitals  V>i  { i = 1,  2, ...,  N }.  The  excited  states  were  determined  by 
single  excitations  within  this  determinant.  The  Fock  operator  F’  for  the  solvated 
molecule  was  defined  by 


where  g has  been  defined  in  equation  (2.27)  and  Fo  is  the  Fock  operator  for  the 
solute  in  the  gas  phase  within  the  INDO  approximation.  Comparing  the  terms 
in  equation  (2.32)  with  equation  (2.25),  the  energy  of  the  excited  state  under 
consideration  has  to  be  corrected  by  the  term  r where 


To  evaluate  this  correction,  the  only  additional  quantity  required  is  the  dipole 
moment  of  the  excited  state  of  interest.  For  the  numerical  calculations,  the 
dielectric  constant  and  refractive  index  were  obtained  from  experiment,  and  the 


2.33 
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dipole  moments  of  the  ground  and  excited  states  of  the  solute  calculated  using 
the  INDO  method.  The  only  other  variable  that  has  to  be  determined  is  the  cavity 
radius  'a',  and  the  methods  used  are  outlined  in  the  following  section. 


In  their  studies  on  the  dissociation  of  HF  in  water  using  the  reaction  field 
method,  Karelson  Katritzky,  and  Zemer  (20)  calculated  the  cavity  radius  from  the 
molar  refraction  according  to 


where  n and  Na  are  the  refractive  index  of  the  solute  and  Avagadro’s  number. 
M and  p are  the  molecular  mass  and  the  mass  density.  However,  this  method 
suffers  from  a major  drawback  in  that  solutes  whose  refractive  indices  are  close 
to  unity  are  calculated  to  have  negligible  radii.  A good  example  is  acetylene;  the 
refractive  index  is  1.00051  and  the  radius,  according  to  equation  (2.34),  is  0.178A. 

Another  method  to  evaluate  radii  is  from  effective  volumes  of  organic  sub- 
stituent groups  proposed  by  Bondi  (21).  The  values  for  'a'  are  realistic,  but  when 
used  in  the  evaluation  of  the  spectral  shifts  as  outlined  in  the  previous  section, 
gave  results  in  poor  agreement  with  experiment. 

Calculation  of  the  cavity  radius  from  the  mass  density  is  another  possibility. 
The  radius  'a'  is  calculated  according  to 


2.5  Determination  of  the  Cavity  Radius  ' a ' 


2.34 


2.35 
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The  spectral  shifts  calculated  using  this  method  are  in  much  better  agreement 
with  experiment  than  those  obtained  from  Bondi’s  radii. 

Fitting  a scaled  cavity  radius  obtained  from  equation  (2.35)  to  the  calculated 
spectral  shifts  indicated  that  a radius  that  is  1.5  times  that  calculated  from  the  mass 
density  gave  the  best  fit  for  ortho-nitroaniline.  It  was  therefore  decided  to  use  this 
technique  to  obtain  the  cavity  radius  for  the  calculation  of  the  solvent  induced 
spectral  shifts.  The  values  used  for  the  cavity  radius  are  given  in  Table  2.1. 


TABLE  2.1  Cavity  radius  (A)  for  the  solute  molecules. 


SOLUTE 

RADIUS  (A) 

o-nitroaniline 

5.08 

m-nitroaniline 

5.54 

p-nitroaniline 

5.51 

p-nitrophenol 

5.50 

pyridine 

4.76 

phenol 

5.01 

2.6  Results  of  Calculations  and  Disussion 

Using  the  method  outlined  in  sections  2.4  and  2.5,  the  shifts  in  the  electronic 
spectra  of  a series  of  molecules  in  various  solvents  were  determined.  The  results 
are  tabulated  below  along  with  the  values  for  the  cavity  radii  (in  A)  used  in  the 
calculations.  The  geometries  of  the  solute  molecules  were  obtained  from  X-ray 
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data,  and  the  dielectric  constant  and  refractive  index  were  taken  from  standard 
tables  (22)  and  are  listed  in  Table  2.2. 


TABLE  2.2  Dielectric  constant  (e)  and  Refractive  index  (n)  for  the  solvents. 


SOLVENT 

e 

n 

n-hexane 

1.890 

1.3751 

n-pentane 

1.844 

1.3575 

n-octane 

1.948 

1.3974 

cyclohexane 

2.023 

1.4266 

CS2 

2.641 

1.6319 

ethanol 

24.3 

1.3661 

methanol 

32.63 

1.3288 

water 

78.0 

1.33 

The  active  space  for  the  singles  only  configuration  interaction  scheme  was 
confined  to  the  five  highest  occupied  and  five  lowest  unoccupied  molecular 
orbitals.  A calculation  using  a full  singles  only  space  for  pyridine  in  hexane 
and  in  ethanol  gave  a spectral  shift  improved  by  only  ten  wavenumbers.  It  was 
thus  decided  to  use  the  previously  defined  active  space  of  ten  molecular  orbitals 
in  all  the  calculations  as  the  results  obtained  from  a full  singles  only  space  did 
not  justify  the  time  and  effort  required  for  their  evaluation. 

The  origins  of  the  calculated  and  experimental  shifts  for  all  the  results  given 
in  this  chapter  are  the  gas  phase  transition  values  obtained  using  the  INDO  scheme 
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and  experiment  respectively.  These  values  are  given  in  Table  2.3.  A negative  sign 
for  the  solvent  induced  shift  indicates  a red  shift  with  respect  to  the  corresponding 
origin. 


TABLE  2.3  Gas  phase  spectra  of  the  solute  molecules. 


SOLUTE 

CALCULATED 

EXPERIMENTAL 

TRANSITION  (cm  -1) 

TRANSITION  (cm  4) 

o - nitroaniline 

27,400 

28,200  (ref  23) 

m - nitroaniline 

28,710 

30,900  (ref  23) 

p -nitroaniline 

31,699 

34,480  (ref  23) 

p - nitrophenol 

33,302 

37,310  (ref  23) 

pyridine 

38,194 

38,350  (ref  24) 

phenol 

36,950 

36,352  (ref  25) 

TABLE  2.4  Shifts  in  the  spectrum  of  o - nitroaniline. 


SOLVENT 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm'1)  (ref  23) 

ethanol 

-2543 

-3321 

n - hexane 

-2040 

-1590 

n - pentane 

-1955 

-1500 

cyclohexane 

-2282 

-1700 

n - octane 

-2148 

-1690 

CS2 

-3217 

-2680 
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TABLE  2.5  Solvent  shifts  in  the  7r— >7r*  spectrum  of  m - nitroaniline. 


SOLVENT 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm1)  ( ref  23) 

n - hexane 

-2554 

-1920 

n - pentane 

-2438 

-1830 

cyclohexane 

-2867 

-2100 

CS2 

-4144 

-2970 

TABLE  2.6  Shifts  in  the  7r— >7r*  spectrum  of  p - nitroaniline. 


SOLVENT 

CALCULTED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm'1)  (ref  23) 

n - hexane 

-3688 

-3130 

n - pentane 

-3529 

-2810 

cyclohexane 

-4136 

-3410 

CS2 

-5914 

-5280 
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TABLE  2.7  Shifts  in  the  tt— »7t*  spectrum  of  p - nitrophenol. 


SOLVENT 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm'1)  (ref  23) 

cyclohexane 

-2799 

-2340 

n - hexane 

-2482 

-2100 

n - heptane 

-2562 

-2340 

The  geometries  used  for  the  nitroanilines  and  p-nitrophenol  were  obtained 
from  references  (26)  and  (27)  respectively. 

In  each  of  the  four  cases  presented  above,  the  calculated  shifts  parallel  the 
experimental  values  for  each  solute  in  the  various  solvents.  With  the  exception 
of  o-nitroaniline  in  ethanol,  the  calculated  shifts  are  overestimated  by  approxi- 
mately 300  to  1000  wavenumbers  towards  the  red  side  of  the  electromagnetic 
spectrum.  This  agreement  with  experiment  is  considered  to  be  excellent  within 
the  approximations  of  the  model. 

The  calculated  red  shift  for  the  solvation  of  o-nitroaniline  in  ethyl  alcohol  is 
lower  than  the  experimental  value.  In  order  to  rationalize  this  observation,  the 
basic  assumptions  of  the  model  have  to  be  reviewed.  The  reaction  field  equations 
were  derived  based  on  the  theory  of  weak  interactions  between  the  solute  and 
the  solvent.  Furthermore,  the  solvent  was  considered  to  be  isotropic.  These 
assumptions  hold  more  closely  when  the  solvent  is  nonpolar,  in  which  case  the 
intermolecular  interactions  are  small  and  the  specific  arrangement  of  the  solvent 
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TABLE  2.8  Shifts  in  the  7r — » n* spectrum  of  pyridine. 


SOLVENT 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm1) 

water 

290 

713  (ref  28  ) 

ethanol 

347 

560  ( ref  29  ) 

n-hexane 

27 

1491  ( ref  30  ) 

cyclohexane 

32 

1491  ( ref  30  ) 

methanol 

313 

1176  (ref  30  ) 

TABLE  2.9  Shifts  in  the  n— >n*  spectrum  of  phenol. 


SOLVENT 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm1) 

SHIFT  (cm1) 

ethanol 

69 

279  (ref  31  ) 

methanol 

75 

278  (ref  32  ) 

water 

84 

685  (ref  33  ) 

cyclohexane 

7 

413  (ref  33  ) 

hexane 

6 

548  (ref  34) 

X-ray  geometries  for  pyridine  and  phenol  were  obtained  from  references  (35) 
and  (36)  respectively. 

The  influence  of  the  solvent  on  the  ground  and  excited  states  of  the  solute 
molecule  is  responsible  for  the  observed  spectral  shift.  If  the  ground  state  of  the 
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solute  is  polar,  it  will  be  stabilised  by  a polar  solvent.  As  the  time  interval  of 
the  electronic  transition  is  small,  the  solvent  cage  does  not  have  sufficient  time 
to  rearrange  itself  to  accommodate  the  changed  dipole  of  the  solute  excited  state. 
If  the  dipole  moment  or  the  largest  Cartesian  component  of  the  dipole  is  larger 
and  parallel  to  the  corresponding  property  of  the  molecule  in  the  ground  state, 
the  excited  state  solute  will  be  stabilized  to  a greater  extent  than  the  ground  state, 
resulting  in  a red  shift  in  the  spectra.  Clearly,  a nonpolar  solvent  in  the  above 
case  will  yield  a blue  shift. 

In  cases  where  the  excited  state  has  a lower  dipole  than  the  ground  state,  the 
solvent  will  not  stabilize  the  excited  state  as  much  as  the  ground  state,  resulting  in 
a blue  shift  in  the  electronic  spectrum.  If  the  solvent  were  nonpolar,  the  direction 
of  the  shift  will  be  in  the  opposite  direction,  thus  causing  a red  shift  in  this  case. 
The  situations  discussed  in  this  and  the  preceding  paragraph  apply  to  tv  — > tv* 
transitions. 

The  n-orbital  is  responsible  for  the  hydrogen  bond  formation  between  a polar 
molecule  and  a polar  solvent.  During  an  n — > tv*  transition,  an  n-electron  is 
removed  from  the  non-bonding  orbital  and  promoted  to  the  n*  orbital  thereby 
destroying  the  hydrogen  bond.  Thus  the  excited  energy  is  raised  by  an  amount 
equal  to  the  strength  of  the  broken  bond  causing  a blue  shift. 

The  spectra  studied  for  both  pyridine  and  phenol  are  of  the  -k  — * tv*  type,  and 
according  to  the  discussion  presented  in  the  preceding  paragraphs  should  have  red 
shifts  in  polar  solvents  and  blue  shifts  in  nonpolar  environments.  Both  molecules 
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are  predicted  to  have  hypsochromic  shifts  in  all  the  solvents  studied  in  agreement 
with  experiment. 

The  polar  molecules,  water,  ethanol  and  methanol  have  acidic  hydrogens 
which  can  protonate  the  nitrogen  atom  in  pyridine.  Thus,  the  solute  that  is 
being  investigated  experimentally  is  no  longer  pyridine,  but  the  pyridinium  ion 
which  is  now  stabilized  by  strong  hydrogen  bonding  with  the  solvent  molecules 
resulting  in  a blue  shift  in  the  7r  — > n*  spectrum.  The  blue  shift  is  predicted  to 
increase  with  increasing  acidity  of  the  solvent,  an  effect  that  has  been  observed 
experimentally(30). 

Recent  multiphoton  resonance  spectral  studies  on  the  stepwise  hydration  of 
phenol  (35)  show  that  the  lowest  n — ► n*  spectrum  of  phenol  / (H20)i  gives  a red 
shift  indicating  that  the  phenol  is  behaving  as  an  acid.  Further  hydration  shows  a 
continuous  blue  shift  in  the  spectrum  of  the  solute  due  to  the  stabilization  of  the 
phenolic  anion  by  the  solvent  molecules  which  are  capable  of  hydrogen  bonding, 
an  effect  very  similar  to  the  case  of  pyridine. 

Thus,  the  general  guidelines  for  the  prediction  of  the  spectral  shifts  in  polar 
and  nonpolar  solvents  presented  earlier  in  this  section  are  helpful  only  if  the 
solvent  does  not  change  the  nature  of  the  solute. 

2.7  Conclusions 

The  continuum  model  works  very  well  for  solutes  in  nonpolar  solvents, 
provided  a judicious  choice  is  made  for  the  cavity  radius.  In  the  case  of  polar 
solvents,  the  predictive  power  of  the  method  is  not  very  useful,  in  that  the  value 
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calculated  for  the  shift  is  not  in  very  good  agreement  with  experiment.  However, 
the  direction  of  the  shift  is  always  reproduced  correctly. 

Three  immediate  improvements  to  the  method,  as  outlined  in  this  chapter  are 
clear.  The  first  would  be  to  determine  the  cavity  radius  using  a more  reasonable 
algorithm  as  the  spectral  shift  is  strongly  dependent  on  this  parameter.  If  the 
value  of  this  radius  is  too  large,  the  solvent  effect  is  underestimated,  sometimes 
resulting  in  shifts  that  are  nearly  zero.  A cavity  which  is  too  small  overestimates 
the  shift  giving  results  that  are  physically  unrealistic.  Examples  of  such  cases  are 
discussed  in  Chapter  4 where  a suitable  choice  for  the  radius  has  to  be  made  in 
order  to  include  previously  determined  solvent  shells  explicitly.  Furthermore,  the 
same  value  for  the  radius  has  been  used  for  the  ground  and  excited  states.  The 
excited  state,  especially  when  the  solvent  is  polar,  would  be  expected  to  have  a 
larger  radius  as  the  solvent  would  stabilize  the  more  diffuse  charge  density  and 
hence  the  larger  dipole.  It  may  thus  be  necessary  to  use  different  values  for  the 
cavity  radii  when  treating  the  ground  and  excited  states. 

In  all  the  applications  presented  in  this  chapter,  the  solute  molecules  have 
permanent  dipoles.  The  reaction  field  theory,  as  presented  here,  is  applicable 
only  to  such  molecules  and  this  restricts  any  investigations  on  nonpolar  molecules 
like  benzene.  However,  such  molecules  can  be  studied  using  this  formulation  if 
the  two  second  order  terms,  namely  Ei(2)  and  the  dipersion  energy  E3®  ignored 
here  are  explicitly  considered.  The  term  Ei(2),  which  is  the  energy  contribution 
due  to  the  polarization  of  the  solute  by  the  solvent,  will  be  the  leading  term  in 
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the  multipolar  expansion  of  the  interaction  energy  for  nonpolar  systems  and  the 
dispersion  energy  term  will  also  be  appreciable.  A second  improvement  would 
then  involve  extension  of  the  method  to  include  these  two  terms. 

The  third  improvement  is  to  explicitly  consider  the  structure  of  the  solvent 
molecules,  particularly  those  in  the  vicinity  of  the  solute  molecule.  After  con- 
structing these  solvent  shells,  the  solvated  solute  complex  can  then  be  immersed 
in  the  continuum  solvent,  the  radius  of  the  cavity  for  the  complex  having  been 
suitably  determined.  The  solvation  shells  would  introduce  solvent  anisotropy,  an 
effect  overlooked  by  the  continuum  fluid  model.  A method  for  the  construction 
of  such  solvation  shells  is  presented  in  the  next  chapter. 


CHAPTER  3 

SEARCH  FOR  THE  GLOBAL  MINIMUM 
ON  A POTENTIAL  ENERGY  SURFACE 


3.1  Introduction 


To  obtain  the  most  stable  structure  of  solvent  molecules  around  a solute,  it  is 
necessary  to  find  the  global  minimum  on  the  potential  energy  surface  (PES)  that 
describes  the  solute-solvent  interactions.  This  search  can  be  accomplished  in  two 
ways,  namely 

1.  Sequential  search  — gradient  based  methods,  and 

2.  Random  search  techniques. 

The  sequential  search  algorithms  require  the  evaluation  of  the  gradients  fi 
at  points  r,  on  the  surface  (38)  which  are  generated  from  an  arbitrarily  chosen 
starting  point  ro  according  to  the  equation 

r,+i  = ri  - a/,#"1  3.1 

where  i = 0,  1,  2,  . . . , N 
a is  the  step  size  on  the  surface 
H is  the  second  derivative  of  the  energy  matrix  and 
N is  the  number  of  steps  allowed  in  the  search. 
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For  a simple  PES  with  one  or  very  few  minima  which  are  sparsely  distributed, 
the  gradient  search  is  efficient  in  locating  these  minima.  However,  for  systems 
with  many  degrees  of  freedom,  like  those  found  in  problems  dealing  with  solva- 
tion, the  gradient  method  fails.  This  is  because  the  search  can  be  trapped  in  any 
one  of  the  local  minima  and  the  algorithm  does  not  have  a mechanism  by  which 
it  can  move  out  of  such  minima. 

Furthermore,  ro  defines  the  point  on  the  the  PES  at  which  the  search  is 
initiated.  The  success  of  the  algorithm  depends  on  the  choice  of  this  point,  and, 
in  the  case  of  surfaces  with  multiple  minima,  the  scheme  is  efficient  only  if  ro 
is  in  the  neighborhood  of  the  global  minimum.  If  this  condition  is  not  satisfied, 
the  local  minimum  closest  to  the  choice  of  ro  will  be  calculated.  However,  the 
advantage  in  using  gradient  type  methods  is  that  they  are  computationally  less 
time  consuming. 

Random  search  techniques,  on  the  other  hand,  do  not  have  these  disadvan- 
tages. The  potential  energy  surface  is  scanned  at  random,  the  new  point  sampled 
being  a function  of  a predetermined  step  size  of  the  various  degrees  of  freedom. 
The  method  of  Simulated  Annealing  (39,40)  as  adapted  by  Donnelly  (41)  is  an 
example  of  a random  search  technique  which  is  applicable  to  problems  in  intra- 
and  intermolecular  interactions.  As  the  systems  studied  in  this  thesis  deal  with 
the  latter  category,  the  outline  of  the  method  presented  here  will  emphasize  only 
this  aspect  of  the  algorithm. 


36 


3.2  The  Simulated  Annealing  Method 


A very  successful  method  by  which  the  equilibrium  properties  of  systems 
can  be  evaluated  is  the  Metropolis  Monte  Carlo  algorithm  (42).  Thermodynamic 
observables  are  a statistical  average  of  the  property  measured,  and  this  average 
is  over  an  order  of  the  Avagadro’s  number  in  the  constituent  particles.  Computa- 
tionally, such  calculations  are  impossible  and  it  is  therefore  necessary  to  simulate 
the  bulk  properties  from  a small,  but  physically  meaningful  number  of  particles. 
The  choice  of  this  number,  say  N,  should  be  so  as  to  minimize  the  fluctuations 
in  the  calculated  value  of  the  property. 

The  N particles  are  assumed  to  interact  with  each  other  and  this  interaction 
is  described  by  an  interparticle  potential  energy  function  U(r,fl).  The  variables  r 
and  Q.  represent,  respectively,  the  distance  separating  the  particles  and  any  other 
coordinates  such  as  the  angles  of  orientation  that  the  potential  energy  might  depend 
on.  The  interparticle  interaction  energy  AE  is  then  given  by 


In  the  Monte  Carlo  method,  the  N particles  are  placed  in  a cube  whose 
volume  V has  been  chosen  to  reproduce  experimental  data  on  the  mass  density  of 
the  system  to  be  simulated.  The  initial  configuration  of  the  particles  is  arbitrary 
and  this  flexibility  is  one  of  the  tremendous  advantages  of  this  method.  Using 
equation  (3.2),  AE  for  the  initial  arrangement  is  evaluated. 


3.2 
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A particle  i at  position  r,  = (xj,yj,Zj)  is  chosen  at  random  and  moved  according 
to 

x\  = z,  + dn\ 
y'i  = Vi  -f  dn2 

z[  = z,  + dna  3.3 


ft;  = ft'  + fnk 

where  d and  f are  previously  chosen  step  sizes  for  the  changes  in  the  displacement 
rand  ft  variables.  The  set  { ni } i=l,2, are  random  numbers  in  the  interval  [-1,1]. 

If  the  move  takes  the  particle  out  of  the  cube,  the  particle  is  brought  back  into 
the  volume  V from  the  direction  which  is  the  mirror  image  of  the  displacement 
that  took  the  particle  out  of  the  cube.  This  requirement  reduces  the  surface  effects 
and  maintains  the  mass  density  at  its  predetermined  value. 

The  interaction  energy  in  the  new  configuration  AE*  is  now  calculated  using 
equation  (3.2).  If  AE  is  less  than  the  value  calculated  for  the  previous  configura- 
tion, the  new  arrangement  of  the  particles  is  accepted  and  the  process  of  choosing 
a particle  and  moving  it  to  a newer  configuration  repeated.  In  cases  where  AE  is 
greater  than  the  previous  value  for  the  interaction  energy,  a probability  p,  given  by 

p = exp  [-  (A E'  - AE)  /kT]  3.4 

is  calculated,  k and  T refer  to  the  Boltzmann  constant  and  the  thermodynamic 
temperature  respectively.  If  p is  greater  than  a random  number  belonging  to  the 
interval  [0,1],  the  move  of  the  particle  to  the  new  configuration  is  accepted  and 
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the  algorithm  proceeds  as  descibed  above.  On  the  other  hand,  if  p is  less  than  the 
random  number,  the  new  configuration  is  rejected  and  the  particle  brought  back 
to  its  position  prior  to  the  move.  A particle  is  again  chosen,  moved  according 
to  equation  (3.3),  and  the  method  repeated  until  no  further  configurations  are 
accepted.  When  this  criterion  is  satisfied,  the  system  is  said  to  have  reached  its 
equilibrium  configuration. 

The  efficiency  with  which  this  minimum  is  reached  according  to  the  Metropo- 
lis algorithm  will  depend  on  the  number  of  moves  allowed  for  the  displacement 
of  the  particles.  The  equilibrium  configuration  of  the  particles  will  be  determined 
more  accurately  if  the  number  of  moves  allowed  is  very  large  in  the  random  walk 
process. 

The  simulated  annealing  method  makes  use  of  the  Monte  Carlo  algorithm 
outlined  above.  The  main  difference  is  that  the  probability  p is  now  calculated 
according  to 

p = exp  [-  (A E'  - A E)  /T*]  3.5 

where  T*  is  now  a parameter  which  has  the  same  units  as  energy.  For  a given 
value  of  T*  the  potential  energy  surface  is  scanned  in  a finite  number  of  moves 
using  the  random  walk  process  described.  T*  is  then  varied  using  an  annealing 
factor  /?  as 

T*+1=PT*,  0</?<l  3.6 

The  index  i refers  to  the  number  of  annealing  steps  allowed.  The  entire  search 
process  is  repeated  with  the  new  value  of  T*.  As  T*  decreases,  areas  of  the  the 
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surface  closer  to  the  minimum  are  scanned  and  if  any  minima  had  been  missed 
by  the  search  using  the  previous  value  of  T*,  the  algorithm  can  now  lock  itself 
onto  a lower  minimum.  As  the  number  of  cycles  allowed  for  the  annealing  step 
increases,  the  search  for  the  global  minimum  becomes  more  efficient  as  minima 
of  lower  energy  can  be  searched.  This  flexibility  in  being  able  to  anneal  the  PES 
is  one  of  the  assets  of  the  simulated  annealing  method. 

Unlike  gradient  search  techniques  in  which  the  displacement  or  steps  are  on 
the  potential  surface,  the  random  displacements  in  annealing  enable  the  search  to 
tunnel  out  of  any  local  minima  in  which  the  algorithm  could  have  got  trapped. 
The  decrease  in  T*,  which  allows  narrower  areas  of  the  surface  to  be  searched 
coupled  with  the  tunneling  advantages  assures  the  evaluation  of  the  global  mini- 
mum provided  the  algorithm  has  been  allowed  a large  enough  number  of  random 
moves  for  each  value  of  T*. 

3.3  Choice  of  Potential  Energy  Function 


In  order  to  apply  the  annealing  method  described  in  the  previous  section,  the 
interaction  energy  has  to  be  evaluated.  This  requires  a prudent  choice  of  U(r,fi). 
A quantum  mechanical  calculation  is  impractical  as  the  energy  has  to  be  evaluated 
for  changes  in  every  intersystem  coordinate.  The  energy  of  interaction  AE  can 
be  expressed  as  a series  according  to  (5) 

ae  = ^u  (r,,  7j,  nufij)  + S^(ri,ri,r*,nilnilnt)  + ... 

i,j,k 


3.7 
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where  U(ri,  . . . ,rm,  Q,\,  . . . ,f)m)  is  the  m-body  contribution  to  the  interaction 
energy.  The  various  terms  in  the  above  equation  refer  to  summations  over  the 
atomic  centers  of  the  constituent  systems. 

Application  of  equation  (3.7)  to  a problem  that  deals  with  molecu- 
lar solvation  on  a macroscopic  scale  is  not  computationally  feasible.  The 

main  reason  is  the  number  of  terms  that  have  to  be  evaluated.  Even 

if  a simulation  of  the  macroscopic  system  is  executed  with  N molecules, 
this  problem  still  remains  as  each  m-body  term  will  have  Nm  calcula- 
tions. This  number  is  astronomical  even  for  small  N for  which  the  se- 
ries expansion  of  AE  has  higher  order  terms  in  the  interaction.  As  such, 
it  is  necessary  to  truncate  the  series  at  a reasonable  order  in  the  inter- 
action which  incorporates  the  physics  of  interest.  Expressing  AE  as  a 
sum  of  interactions  over  pairs  of  atomic  centers  is  a method  commonly 
used. 

If  the  distance  Rab  separating  the  systems  A and  B is  much  larger  than 
the  dimensions  of  the  interacting  species  itself,  the  energy  of  interaction  can 
be  evaluated  using  the  theory  of  weak  interactions.  Within  this  theory,  AE 

is  expressed  as  a multipolar  expansion  in  the  inverse  powers  of  Rab  as 


(5) 
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3.8 


n= 1 


The  coefficients  Cn  contain  the  physics  of  the  interaction.  For  example  when 
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n=l,  Ci  is  given  by 


Ci  = qAVB 


3.9 


which  is  the  coefficient  of  the  coulombic  term.  qA  and  qe  are  the  charges 
on  systems  A and  B respectively.  The  higher  order  coefficients  contain 
dipole  moments,  polarisabilities  and  higher  order  moments  of  the  interacting 
species. 

If  A and  B are  molecular  systems  composed  of  a and  b atomic 
centers  respectively,  the  pairwise  interaction  approximation  expresses  AE 
as 


Clearly,  the  infinite  sum  in  the  above  equation  is  a restriction  to  the  eval- 
uation of  the  interaction  energy.  It  is  thus  necessary  to  truncate  equa- 
tion (3.10)  to  a finite  number  of  terms.  Many  successful  applications  us- 
ing the  series  of  approximations  outlined  above  are  available  in  the  literature 
(43,44). 

The  coefficients  Cn  for  truncated  pairwise  potentials  can  be  obtained 
in  many  ways.  One  method  is  to  use  the  electronic  moments  obtained 
from  experiment  together  with  classical  or  quantum  mechanical  charge  densi- 
ties. Such  a method  suffers  from  a serious  drawback  in  that  the  static  mo- 
ments used  do  not  always  simulate  the  interactions  accurately  and  the  ex- 
perimental values  contain  many  body  effects  which  have  been  neglected  in 
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3.10 


n=l  i = l j=l 
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the  pairwise  approximation.  For  example,  the  description  of  the  interac- 
tion of  a positive  ion  with  a neutral  molecule  should  be  able  to  account  for 
the  change  in  the  polarizabilities  of  the  interacting  systems  due  to  charge 
transfer  effects,  a phenomenon  that  is  neglected  when  static  moments  are 
used. 

Another  method  is  to  fit  a multipolar  expansion  to  thermodynamic  data 
for  the  systems  studied.  A good  example  of  such  a parameterisation  is 
the  potential  energy  function  determined  for  water  by  Jorgensen  et  al  (45). 
The  Bernal  and  Fowler  model  (46)  is  used  to  describe  the  water  mole- 
cule, according  to  which  a molecule  of  water  is  comprised  of  four  copla- 
nar  points,  namely  the  three  atoms  and  a pseudocenter  M on  the  bisector 
of  the  HOH  angle.  This  point  is  at  a distance  of  0.15A  from  the  oxygen 
atom. 

The  interaction  energy  AE  is  given  by 


The  summations  are  over  the  centers  on  the  water  molecule,  q*  is  the  quantum  me- 
chanical charge  density  on  center  i and  Ay  and  By  are  coefficients  determined  as 
described  below.  Using  equation  (3.11),  a Monte  Carlo  simulation  of  the  thermo- 
dynamic properties  of  water  like  the  specific  heat  and  isothermal  compressibility 
is  carried  out.  Ay  and  By  are  varied  to  reproduce  these  quantities  as  closely  as 
possible  to  experimental  data.  The  values  obtained  for  the  parameters  used  in 
equation  (3.11)  are  given  in  Table  3.1. 


3.11 
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TABLE  3.1  Parameters  for  equation  (3.11) 


PARAMETER 

VALUE 

r(OH  A) 

0.9572 

r (OM  A) 

0.15 

Angle  HOH,  deg 

104.52 

A ( kcal  A6/mol) 

610.0 

B x 1 O'3, kcal  A12/mol 

600.0 

q (0) 

0.0 

q (H) 

0.52 

q(M) 

-1.04 

The  use  of  such  potential  energy  functions  is  costly  as  simulations  are  required 
for  every  system  that  is  to  be  studied.  The  transferability  of  these  parameters  for 
the  study  of  mixed  systems,  that  is,  systems  with  more  than  one  chemical  species, 
is  questionable  as  the  specificities  of  the  interactions  in  mixed  systems  have  to  be 
considered.  Furthermore,  a set  of  parameters  for  each  pair  of  interacting  atoms  is 
required  to  define  the  potential  for  a variety  of  systems. 

A very  practical  approach  to  parameterising  a potential  energy  function  that 
can  be  applied  to  a variety  of  systems  was  presented  by  Fraga  (47).  Instead  of 
determining  parameters  for  every  pairwise  interaction,  the  atoms  in  the  interacting 
systems  are  assigned  to  classes  based  on  the  average  one  center  contribution  to 
the  Hartree-Fock  energy  called  the  Molecular  Orbital  Valency  State  (MOVS)  (48). 


44 


The  energy  of  interaction  AE  is  calculated  as 

= 5 ^{qaqb/Rab  - \ [faOiaql  + fb<*bq2a]  /Kb 

a,b 


3 

2 


( fafbdtaOib / ( fa<Xa/na + (/t<Wnfc)1/2 


3.12 


where  the  summation  indices  (a,b)  are  now  over  the  classes  to  which  the  atoms 
belong.  qx,  ax,  nx,  fx  and  gx  (x=a,b)  are  respectively,  the  electronic  charge,  the 
static  polarisability,  the  effective  number  of  electrons  and  parameters  for  class  x. 
qx  is  obtained  from  quantum  mechanical  calculations  (49)  and  the  polarisability  is 
calculated  by  extrapolating  Hartree-Fock  data  for  the  atom  and  its  ions  (50)  that 
have  been  assigned  to  class  x.  nx,  the  effective  charge  is  given  by 

nx  — Zx  — qx  3.13 

where  Zx  is  the  atomic  number  of  the  atom  that  has  been  assigned  to  class  x. 
fx  and  gx  are  obtained  by  fitting  equation  (3.12)  to  ab-initio  results  of  Clementi 
et  al  (49)  for  the  interaction  of  twenty  six  amino  acids  with  water.  Table  3.2 
gives  the  parameters  for  the  atoms  (47)  and  ions  (51)  and  a description  of  their 
environments. 


TABLE  3.2  Parameters  for  Atoms  and  Ions. 
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ENVIRON- 

MENT 

- SH 

-CH.CH2  ; 

-CH3 

- CH, 

AROMATIC 

C 

- NH2 

- COOH  | 

-CH3 

- CH2, 
ALPHA  TO 
S 

-CH2  j 

AROMATIC 

C 

< ~ 
ci  £ 

li 

bO 

69.3 

162 

22.3 

131 

1.28 

5.53 

536 

21.2 

508 

948 

.85 

ro 

0.1 

0.1 

0.35 

0.95 

oro 

© 

© 

oro 

oro 

< 

'o 

0.42 

0.36 

0.36 

0.34 

0.33 

0.29 

3.60 

3.20 

2.75 

2.20 

cr 

0900 

0.204 

0.205 

0.253 

0.266 

0.404 

-0.608 

ro 

O 

2 

9 

ro 

00 

ro 

9 

oo 

rs 

9 

CLASS 

23 

i 

<s 

vO 

- 

T}- 

VO 

r~~ 

r-~ 

2 

X 

u 

5 

TABLE  3.2  (continued) 
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TABLE  3.2  (continued) 
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ENVIRON- 

MENT 

ADJACENT 
TO  NH2 

AROMATIC, 

CO 

CARBO- 

XYLIC 

AMIDE 

AMINO 

AROMATIC, 
WITH  H 

AROMATIC, 

WITHOUT 

H 

X 

Vi 

< ~ 
t*  s 

73  is 
u 5 

* E 

bo 

982 

385 

1 

667 

1111 

245 

998 

362 

j 1035 

2.60 

oro 

oro 

oro 

0.15 

0.95 

oro 

| 0.35 

< 

C 

1.30 

1.20 

Oil 

2.35 

2.10 

1.90 

1.50 

3.15 

CT 

0.310 

0.423 

uso 

-0.630 

-0.554 

-0.473 

-0.317 

0.123 

CLASS 

24 

26 

<N 

22 

ATOM 

U 

z 

Vi 

TABLE  3.2  (continued) 
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ENVIRON- 

MENT 

.E  g 

1 8 
O 0 

CO,  IN 
COOH 

CO, 

ATTACHED 
TO  RING 

0 

fM 

1 

0 

C4 

1 

v O j 

< ~ 
C £ 

T3  % 

Si 

60 

555 

476 

268 

548 

547 

I 563 

t 1415 

Um  | 

0.21 

oro 

080 

0.95 

060 

1 5.1 

001  j 

< 

e 

1.45 

1.25 

1.20 

1.65 

0.32 

| 0.15 

| 0.95 

o* 

-0.539 

-0.409 

o 

op 

in 

o 

i 

-0.622 

0.311 

g 

+ 

8 

+ 

00 

00 

hs 

o 

<s 

PO 

3 

ON 

CN 

CO 

<N 

CO 

CO 

u 

2 

+ 

p 

O 

o 

X 

cd 

z 

* 

§ 
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Both  potential  functions  mentioned  above  were  used  to  evaluate  the  interaction 
energy  of  simple  systems.  The  water  dimer  was  the  first  example  studied  using 
the  simulated  annealing  method.  The  initial  orientation  for  the  system  is  given 
in  Figure  3.1.  The  geometry  used  for  the  water  molecule  was  obtained  from 
reference  (52)  and  the  results  obtained  are  listed  in  Table  3.3.  6 is  the  angle  made 
by  the  OH  bond  of  the  second  molecule  with  the  bisector  of  the  HOH  angle  of 
the  first  molecule. 


B 


Figure  3.1  Initial  (A)  and  final  geometries  (B)  for  the  water  dimer. 
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3.4  Results  and  Discussion  of  the  Simulated  Annealing  Method 


The  starting  configuration  for  the  water  dimer  (Figure  3.1)  used  in  the  previous 
section  is  one  that  is  very  close  to  the  equilibrium  configuration.  The  potential  of 
the  annealing  method  lies  in  the  fact  that  the  results  are  independent  of  the  initial 
arrangement  of  the  molecules.  In  order  to  illustrate  this  aspect  of  the  algorithm, 
the  calculation  on  the  water  dimer  was  repeated  with  the  two  molecules  arranged 
as  shown  in  Figure  3.2. 


Figure  3.2  Initial  random  geometry  for  the  water  dimer. 

The  configuration  shown  above  has  a very  high  energy  as  the  lone  pairs  on 
the  oxygen  atoms  repel  each  other.  A gradient  search  technique  fails  to  give 
a minimum  even  after  many  hours  of  computation  as  the  two  molecules  repel 
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TABLE  3.3.  Calculated  and  experimental  results  for  the  water  dimer 
using  the  Jorgensen  (9)  and  Fraga  (11)  potentials. 


PROPERTY 

EXPERIMENT 

JORGENSEN 

FRAGA 

(52-54) 

POTENTIAL 

POTENTIAL 

r O-H  (A) 

0.9572 

HOH 

104.52 

AE  (kcals/mole) 

-5.44 

-6.28 

-5.07 

e deg 

50-70 

73.9 

66.2 

r O1O2 

2.98 

2.91 

2.89 

The  interaction  energy  AE  we  obtained  using  the  first  potential  function  is 
lower  than  the  value  of  - 6.24  kcals/mole  reported  by  Jorgensen  et  al  (45)  from 
their  Monte  Carlo  simulation,  indicating  that  the  minimum  they  had  evaluated  with 
their  potential  is  not  the  global  minimum  for  the  system.  Further,  the  agreement 
with  experiment  for  the  other  parameters  is  not  nearly  as  good  as  that  obtained 
using  the  Fraga  potential. 

One  area  of  study  presented  in  this  dissertation  is  the  interaction  of  molecules 
with  a wide  range  of  solvent  systems.  It  is  therefore  important  to  have  one  po- 
tential energy  function  which  has  been  consistently  parameterised  and  applicable 
to  as  many  systems  as  possible.  The  comparison  of  results  obtained  for  different 
systems  from  simulations  that  have  been  parameterised  by  various  means  loses 
meaning  in  their  interpretation.  Taking  these  points  into  consideration,  it  was 
decided  to  use  the  Fraga  potential  in  all  the  simulations  to  be  presented  in  this 
thesis.  The  results  of  these  simulations  are  presented  in  the  following  section. 


52 


each  other  very  slowly.  The  simulated  annealing  method  is  very  successful  as 
the  results  given  below  illustrate. 


TABLE  3.4  Results  for  the  water  dimer 


PROPERTY 

GEOMETRY 

GEOMETRY 

EXPERIMENT 

(FIGURE  3.1) 

(FIGURE  3.2) 

(ref  52-54) 

AE  (kcals  /mole) 

-5.07 

-5.07 

-5.44 

r 0i02  (A) 

2.89 

2.89 

2.98 

0 deg 

66 

67 

50-70 

The  calculated  properties  for  the  water  dimer  shown  in  coloumns  2 and  3 
of  table  3.4  show  that  the  results  are  independent  of  the  starting  configuration 
and  in  excellent  agreement  with  experiment.  Thus,  the  initial  configuration  of 
the  systems  to  be  simulated  need  not  be  chosen  beforehand  to  be  close  to  the 
equilibrium  arrangement  of  the  molecules  — a bias  usually  introduced  in  gradient 
search  methods,  (but  not  for  this  case). 

The  water  trimer  was  studied  next.  The  global  minimum  evaluated  corre- 
sponds to  a cyclic  structure  in  which  the  three  oxygen  atoms  form  an  isoceles 
triangle.  Infra-red  studies  (55)  show  the  trimer  to  be  cyclic  with  the  oxygen  atoms 
at  the  vertices  of  an  equilateral  triangle.  A possible  reason  for  the  discrepancy 
between  the  calculated  and  experimental  results  lies  in  the  manner  in  which  the 
interaction  energy  is  evaluated.  The  evaluation  of  the  interaction  energy  as  a sum 
of  pair-wise  terms  neglects  three  body  and  higher  effects,  contributions  which 
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are  present  in  experimental  studies.  Although  these  terms  may  be  comparatively 
small,  they  could  accumulate  to  favour  a different  structure  from  that  obtained 
from  a pair-wise  expansion.  Calculations  on  the  water  tetramer  gives  a cyclic 
structure  as  the  most  stable  configuration  in  agreement  with  other  calculations 
reported  (56-58). 

Another  example  of  a system  stabilised  by  hydrogen  bonding  is  the  acetic  acid 
dimer.  A series  of  minima  were  obtained  and  these  values  are  given  in  Table  3.5. 
The  global  minimum  corresponds  to  an  interaction  energy  of  - 9.46  kcals/mole 
and  the  structure  of  the  dimer  is  given  in  Figure  3.3. 

At  the  global  minimum,  the  carboxylic  groups  in  the  monomers  are  coplanar 
facilitating  maximum  hydrogen  bonding  between  the  carbonyl  oxygen  of  one 
molecule  and  the  hydrogen  atom  of  the  — OH  group  of  the  other.  In  each  of 
the  structures  that  correspond  to  the  local  minima,  the  two  — COOH  groups  are 
not  coplanar  resulting  in  structures  of  higher  energies  due  to  reduced  hydrogen 
bonding. 


TABLE  3.5  Minima  for  the  acetic  acid  dimer 


MINIMUM 

AE  (kcals/mole) 

Local 

-6.87 

Local 

-7.81 

Local 

-8.13 

Local 

-8.29 

Local 

-9.22 

Global 

-9.46 

54 


Figure  3.3  Orientation  of  the  monomers  of  acetic  acid  at  their  global  minimum. 

In  their  studies  on  the  benzene-water  system,  Karlstrom  et  al  (59)  report  two 
minima  — the  global  minimum  being  1.09  kcals/mole  below  the  local  minimum. 
The  annealing  method  gives  very  similar  results  which  are  tabulated  below.  The 
two  sets  of  results  in  Table  3.6  are  not  directly  comparable  as  the  potential  energy 
functions  used  in  their  evaluation  are  different  — the  emphasis  is  on  the  trends 
and  similarities  between  them. 
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TABLE  3.6  Results  for  the  benzene-water  system 


PROPERTY 

SIMULATED 

KARLS TROM  et  al  (ref 

ANNEALING 

59) 

AE  (kcals/mole)  (local 

-1.77 

-1.90 

minimum) 

AE  (kcals/mole)  (global 

-2.66 

-2.99 

minimum) 

0 deg 

0.09 

0.1 

R(A) 

4.75 

3.5 

The  global  minimum  for  this  system  is  shown  in  Figure  3.4.  Both  methods 
give  the  same  orientation  for  the  subsystems,  although  the  intermolecular  distance 
R,  measured  from  the  oxygen  atom  to  the  center  of  the  benzene  ring,  is  different. 


Figure  3.4  Global  minimum  for  the  benzene  — water  system. 
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Clusters  of  molecules  can  be  formed  by  the  stepwise  addition  of  monomers. 
The  simulated  annealing  method  was  applied  to  the  formation  of  such  clusters  or 
microdroplets.  The  global  minimum  for  the  stepwise  formation  of  each  cluster 
was  calculated  and  at  each  step  in  the  buildup  procedure,  all  the  molecules  were 
allowed  to  move  according  to  the  algorithm. 

Based  on  the  solvophobic  theory  of  Sinanoglu  (60),  it  is  possible  to  determine 
the  number  of  molecules  in  a microdroplet,  which  is  defined  as  the  smallest  pos- 
sible cluster  of  the  constituent  molecules  that  approaches  bulk  behaviour.  Using 
the  experimental  values  for  the  surface  tension  and  isothermal  compressibility, 
the  number  of  molecules  in  a cluster  have  been  calculated.  For  water,  a cluster 
of  fifteen  molecules  was  found  to  have  95%  of  the  bulk  surface  tension,  while 
for  benzene  six  molecules  satisfied  this  criterion. 

The  water  and  benzene  systems  were  studied  using  the  annealing  method.  The 
interaction  energy  per  molecule  in  the  cluster,  AEi,  was  evaluated  for  droplets 
of  increasing  size.  The  results  for  benzene  are  tabulated  below  and  the  structure 
for  the  dimer  is  shown  in  Figure  3.5.  The  dimerisation  energy  for  benzene  was 
calculated  to  be  1.50  kcals/mole,  well  within  the  experimental  range  of  1.15  to 
2.07  kcals/mole  (61). 
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TABLE  3.7  Simulation  of  the  microdroplet  of  benzene 


NUMBER  OF  MOLECULES  N 

AE]  (kcals/mole) 

2 

-0.75 

3 

-1.08 

4 

-1.33 

5 

-1.55 

6 

-1.43 

. 

• 

10 

-1.44 

Figure  3.5  Structure  of  the  benzene  dimer  at  its  global  minimum. 
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The  T-shaped  structure  obtained  for  the  dimer  of  benzene  is  in  agreement 
with  experimental  results  (61)  and  quantum  chemical  studies  (62,63).  A plot  of 
AEi  versus  the  number  of  molecules  N is  shown  in  Figure  3.6.  From  the  graph, 
the  number  of  molecules  of  benzene  that  can  form  a cluster  with  no  appreciable 
further  stabilisation  is  five  and  the  structure  of  this  pentamer  is  shown  in  Figure  3.7 


Figure  3.6  Interaction  energy  per  molecule  versus  N for  (benzene)N- 
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Figure  3.7  Structure  of  the  benzene  pentamer. 
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The  curve  shown  if  Figure  3.6  goes  up  for  the  hexamer  and  larger  clusters. 
The  difference  in  the  energy  between  the  minimum  for  the  pentamer  and  that  for 
the  larger  clusters  is  only  0.11  kcals/mole.  This  is  attributed  to  be  an  artifact 
introduced  by  the  potential  energy  function.  The  structure  of  the  pentamer 
(Figure  3.7)  shows  the  molecules  to  be  arranged  perpendicular  to  each  other, 
corresponding  to  a ligand-like  structure,  in  agreement  with  other  theoretical 
calculations  (62). 

A similar  treatment  of  the  water  cluster  yielded  the  results  shown  in  Table  3.8 
and  the  plot  of  the  interaction  energy  per  water  molecule  versus  N,  the  number 
of  molecules,  is  given  in  Figure  3.8 

From  this  figure,  the  stabilisation  of  a cluster  of  water  molecules  shows  no 
further  gain  beyond  N = 16.  The  microdroplet  of  water  is  thus  found  to  have 
sixteen  water  molecules,  in  excellent  agreement  with  the  results  of  Sinanoglu  (60) 
who  reports  a value  of  15.  The  structure  of  this  droplet  is  shown  in  Figure  3.9. 
The  stabilisation  observed  at  N = 5 corresponds  to  the  water  pentamer.  Recent 
theoretical  studies  (64)  show  that  the  first  solvation  shell  of  a water  molecule 
has  four  other  molecules  arranged  tetrahedrally  around  the  central  molecule.  The 
structure  of  the  pentamer  calculated  by  the  annealing  method  is  shown  in  Figure 
3.10. 
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TABLE  3.8  Results  of  the  simulation  of  the  water  cluster 


NUMBER  OF  MOLECULES  N 

AEj  (kcals/mole) 

2 

-2.25 

3 

-4.19 

4 

-5.26 

5 

-5.31 

6 

-5.62 

8 

-6.47 

10 

-6.90 

12 

-7.25 

14 

-7.38 

16 

-8.01 

18 

-7.97 

20 

-7.35 

• 

• 

25 

-7.30 
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Figure  3.8  Interaction  energy  per  molecule  versus  N for  (H20)n. 
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Figure  3.9  Structure  of  the  microdroplet  of  water. 
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Figure  3.10  Structure  of  the  water  pentamer. 

The  interaction  of  positive  ions  with  water  is  of  importance  in  biochemistry. 
Many  theoretical  (65,66)  and  experimental  studies  (67)  on  the  energetics  of  ionic 
solvation,  particularly  for  the  interaction  of  Na+  and  K+  with  water,  are  available 
in  the  literature.  In  the  experimental  investigations,  the  enthalpy  for  the  stepwise 
addition  of  (I^OJn,  n=  1,  . . . , 6,  to  the  sodium  and  potassium  ions  was 
measured.  Cieplak  et  al  (65)  evaluate  the  interaction  energies  for  these  processes 
using  a Monte  Carlo  technique.  The  potential  function  used  in  their  studies 
include,  in  addition  to  the  coulombic,  R~  6 and  R“ 12  contributions,  a term  which 
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depends  explicitly  on  the  transfer  of  charge  from  the  water  molecules  to  the  ion. 
Furthermore,  a three  body  contribution  to  the  interaction  energy  is  also  calculated. 

Parameters  for  ionic  solvation  for  the  Fraga  potential  are  available  (51).  The 
systematic  hydration  energy  per  water  molecule  for  the  stepwise  solvation  of  the 
sodium  and  potassium  ions  were  calculated  using  the  annealing  method.  These 
results  are  given  in  Tables  3.9  and  3.10  (column  2)  along  with  values  obtained 
by  Cieplak  et  al  (65)  in  column  3.  The  corresponding  experimental  values  which 
are  given  in  column  4 have  been  obtained  from  the  measured  enthalpies  corrected 
for  the  kinetic  energy  contribution. 

TABLE  3.9  Interaction  energies  per  water  molecule  (kcals/mole) 
for  the  sodium  ion  with  (water)n 


n 

AEi  (kcals/mole) 

AEi  (kcals/mole) 
(ref  65) 

EXPERIMENT 
(kcals/mole)  (ref 
67) 

1 

-23.89 

-23.01 

-23.40 

2 

-25.30 

-22.35 

-21.30 

3 

-26.04 

-20.43 

-19.27 

4 

-26.69 

-18.53 

-17.93 

5 

-27.05 

-17.24 

-16.54 

6 

-27.39 

-16.33 

-15.47 
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TABLE  3.10  Interaction  energies  (kcals/mole)  per  water  molecule 
for  the  potassium  ion  with  (water)n 


n 

AEj  (kcals/mole) 

AEi  (kcals/mole) 
(ref  65) 

EXPERIMENT 
(kcals/mole)  (ref 
67) 

1 

-18.66 

-17.80 

-17.30 

2 

-20.07 

-16.75 

-16.41 

3 

-20.72 

-15.70 

-15.14 

4 

-21.36 

-14.70 

-14.16 

5 

-21.64 

-13.84 

-13.35 

6 

-21.87 

-13.27 

-12.69 

The  results  in  column  2 of  Tables  3.9  and  3.10  clearly  show  the  inadequacy 
of  the  Fraga  potential  energy  function  and  its  parameters  in  studying  problems 
dealing  with  ionic  solvation.  Not  only  are  the  values  overestimated,  the  trend  in 
the  values  is  completely  opposite  to  that  obtained  from  experiment. 

The  interaction  energy  for  Na+  and  K+  with  one  water  molecule  is  in  excellent 
agreement  with  the  Monte  Carlo  value  and  experimental  results.  The  positive 
charge  on  the  ion  is  not  going  to  remain  localised  on  the  ion  after  the  first  water 
molecule  has  coordinated;  some  of  this  charge  will  be  reduced  by  the  transfer  of 
electron  density  from  the  ligand.  Thus,  the  second  water  molecule  is  no  longer 
interacting  with  a unipositive  ion,  but  an  ion  with  a reduced  charge.  The  multipolar 
expansion,  in  the  weak  interaction  limit,  is  not  able  to  account  for  such  charge 
transfer  terms. 
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For  the  interaction  of  charged  species,  the  coulombic  term  is  certainly  the 
dominant  contribution  to  AE.  As  the  reduction  of  the  charge  on  the  ion  is 
neglected,  the  coulombic  term  is  systematically  overestimated.  Although  the 
short  range  repulsion  term  increases  due  to  the  increase  in  the  number  of  water 
molecules,  it  cannot  adequately  compensate  for  the  large  negative  coulombic 
contribution. 


3.5  Discussion 


The  types  of  systems  to  which  the  simulated  annealing  method  has  been 
applied  in  this  chapter  are  varied.  The  study  includes  clusters  of  polar  and  non 
polar  molecules  and  ion-dipole  interactions. 

The  success  of  the  method  depends  on  the  quality  of  the  potential  energy 
function  used.  If  more  terms  which  describe  the  the  nature  of  the  interactions  are 
included,  the  predictive  power  of  the  method  is  enhanced. 

Within  the  polarisation  approximation,  the  multipolar  expansion  of  the 
interaction  energy  is  exact  only  if  the  summation  is  to  infinite  order.  However, 
it  is  necessary  to  truncate  the  series  to  a finite  numer  of  terms  as  explained  in 
Chapter  1.  In  the  Fraga  potential  function,  the  four  terms  used  to  calculate  AE  are 

1.  the  coulombic  or  long  range  term  that  is  a function  of  R_  ], 

2.  a medium  range  term  which  decays  according  to  R_  4, 

3.  the  short  range  dispersion  term  which  behaves  as  the  inverse  sixth 

power  in  R,  and 
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4.  the  very  short  range  repulsion  term  which  is  a function  of  R 12. 

The  R_  4 term  takes  into  account  the  induction  in  a molecule  caused  by  the 
field  generated  by  the  other  molecules  in  the  vicinity.  Such  a term  is  the  leading 
term  in  the  interaction  of  non-polar  systems  like  benzene.The  success  of  the 
method  in  obtaining  the  T-shaped  structure  for  the  benzene  dimer  is  due  to  the 
inclusion  of  this  term  in  the  evaluation  of  the  interaction  energy. 

For  neutral,  non-polar  systems,  the  leading  term  that  contributes  to  stabil- 
isation is  the  dispersion  term  which  depends  on  R~ 6 as  its  leading  term  and 
is  difficult  to  obtain  from  calculations.  Parameterisation  of  this  term  using  the 
Slater-Kirkwood  approximation  (68)  includes  some  of  the  dispersion  effects. 

The  multipolar  expansion  is  strictly  valid  only  in  the  limit  of  large  intersystem 
separation,  R.  Thus  effects  like  charge  transfer,  which  are  short  ranged  phenom- 
ena, are  not  included.  The  predictive  nature  of  potential  functions  based  on  this 
approximation  for  systems  in  which  such  effects  are  dominant  is  certainly  ques- 
tionable.. The  results  obtained  for  ionic  hydration  (Tables  3.9  and  3.10)  illustrate 
the  limitation  of  the  Fraga  potential.  A careful  study  of  the  results  show  a need 
for  the  incorporation  of  the  charge  transfer  and  other  short-range  terms  into  the 
potential  function  before  it  can  be  used  to  examine  ionic  hydration. 

The  pairwise  summation  of  the  interaction  fails  to  take  three  body  and  higher 
order  effects  into  consideration.  In  the  case  of  dipolar  systems  like  water, 
such  terms  are  important  to  account  for  the  cooperativity  and  the  clustering  of 
the  constituent  molecules.  Simulations  of  other  thermodynamic  properties  like 


69 


free-energy  changes  and  specific  heats,  especially  for  polar  molecules,  have  been 
shown  (69,70)  to  require  the  inclusion  of  such  terms,  in  order  to  obtain  reasonable 
agreement  with  experimental  observations. 

Application  of  the  Fraga  potential  incorporated  into  the  simulated  annealing 
method  has  been  particularly  successful  for  large  systems,  as  illustrated  in  the 
cluster  formation  of  water  and  benzene.  The  microdroplet  sizes  predicted  are  in 
remarkable  agreement  with  those  obtained  from  the  solvophobic  theory. 

The  application  of  the  simulated  annealing  method  using  the  potential  energy 
function  given  in  equation  (3.12)  is  successful  in  finding  the  global  minimum  of 
a molecular  potential  energy  surface.  Its  sucessful  application  to  a wide  range 
of  electronically  neutral  systems  makes  it  particularly  useful  in  finding  the  global 
minimum  on  the  potential  energy  surface  of  solute- solvent  interactions,  thereby 
determining  the  structure  of  solvent  molecules  in  the  solvation  shells  surrounding 
the  solute. 


CHAPTER  4 

ELECTRONIC  SPECTRA  USING  THE 
DISCRETE-CONTINUUM  SOLVENT  MODEL 

4.1  Introduction 

The  anisotropy  of  the  solvent  can  be  taken  into  account  by  assigning  a 
structure  to  the  molecules  that  form  the  solvent  shells  immediately  surrounding 
the  solute.  These  shells  are  referred  to  as  the  discrete  structure  of  the  solvent  and 
the  simulated  annealing  method  described  in  the  previous  chapter  is  well  suited 
in  accomplishing  this  task.  After  the  discrete  solvent  structure  has  been  obtained, 
the  energy  of  the  solute- solvent  cluster  can  be  calculated  using  the  supermolecule 
approach  or  the  point  charge  approximation  to  the  solvent.  The  first  method  has 
been  outlined  in  chapter  one.  In  the  point  charge  model,  the  atoms  in  the  solvent 
molecules  can  be  approximated  as  point  charges,  having  the  same  values  as  those 
used  in  the  simulation.  The  effect  of  these  charges  can  then  be  incorporated 
into  the  solute  Hamiltonian  Hs.  The  derivation  of  the  solute-discrete  solvent 
Hamiltonian,  Hs  _ d>  within  the  point  charge  approximation  for  the  simulated 
solvent  structure  is  given  in  section  4.2. 

Multiphoton  spectroscopy  (MPS)  (71)  is  an  analytical  method  in  which  the 
spectrum  of  a molecule  surrounded  by  a few  solvent  molecules  can  be  obtained. 
The  number  of  solvent  molecules  in  such  clusters  is  usually  in  the  range  of  one 
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to  six.  Using  the  method  presented  in  section  4.2,  the  electronic  spectra  of  a 
series  of  such  clusters  have  been  calculated  within  the  INDO  approximation.  The 
results,  along  with  the  corresponding  MPS  values  are  presented  in  section  4.3. 

To  describe  the  solvent  completely,  each  discrete  solvent  cluster  obtained 
from  the  simulated  annealing  method  is  centered  in  a cavity  having  a suitable 
radius,  which  is  then  embedded  in  the  continuum  solvent.  The  method  used  in 
the  selection  of  a value  for  this  radius  is  presented  in  section  4.4.  The  criterion 
for  building  the  discrete  solvent  cluster  is  that  one  solvent  molecule  at  a time 
is  allowed  to  interact  with  the  solute  molecule.  As  more  solvent  molecules  are 
added,  the  interaction  energy  per  solvent  molecule  AEi  converges  to  an  average 
value.  When  the  decrease  in  AEi  is  less  than  half  a kilocalorie  per  mole,  the 
cluster  buildup  is  considered  to  be  complete  and  the  solvent  molecules  that  form 
such  a cluster  is  said  to  comprise  the  discrete  solvent  structure. 

Treating  the  solute-solvent  cluster  as  a modified  solute,  the  continuum  solvent 
effect  is  incorporated  into  Hs  _ d in  the  same  manner  as  that  presented  in  Chapter 
2.  This  model  is  referred  to  as  the  discrete-continuum  model  and  the  results 
obtained  using  this  model  are  tabulated  in  section  4.4.  Conclusions  based  on  the 
results  presented  in  this  chapter  are  given  in  the  final  section. 

4.2  Derivation  of  the  Hamiltonian  for  the  Solute-Discrete 
Solvent  System 

Let  the  solute  Hamiltonian  be  denoted  Hs.  Suppose  the  M atoms  that  comprise 
the  N solvent  molecules  are  replaced  by  point  charges  { qi  } at  positions  { r*  }, 
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i = 1,  2,  . . . , M.  The  charges  q*  are  the  same  values  as  those  used  in  the 
simulated  annealing  method  to  obtain  the  structure  of  the  solvent  cluster. 

H$_d>  the  Hamiltonian  for  the  interaction  of  the  solute  with  the  discrete 
solvent  structure  is  given  by 


M 


M 


Z A .=  1 r,A  .=1  jGS  r'J 

The  first  term  on  the  right  hand  side  of  the  above  equation  is  the  attraction 


between  the  point  charge  q;  and  the  nuclear  cores  of  the  solute  and  the  second 
term  is  the  electron  - electron  interaction  between  the  solute  electrons  j and  the 


point  charge.  The  summation  index  A represents  the  nuclei  in  the  solute,  Za  is 
the  core  charge  on  center  A and  j runs  over  the  electrons  in  the  solute  molecule  S. 
The  Hamiltonian  H*  for  the  composite  solute-solvent  shell  system  is  given  by 


M 


H-  = Hs  + -Y,E 


ags  i=i 


Za<U 

rAi 


M 


M 


- EEf+iE 

,=i  jes  1 i=i, j=i 


M] 

rXJ 


4.2 


The  energy  E*  is  then  given  by 


E*  =<  *\H*\V  > 


where  'I'  is  the  normalised  electronic  wavefunction  for  the  solute  system.  As  q* 
{ i = 1,  2,  . . . , M } does  not  depend  on 

F = e«  + 5EE-  + 5 E — 

2 1 r Ai  2 r,, 

A£S  .=1  At  »=1, j=l  X] 


4.3 
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where  Enuc.s  is  the  core  repulsion  energy  of  the  solute.  From  the  last  term  of 
equation  4.3,  the  effect  of  the  point  charges  is  to  modify  the  one  electron  part  of 
the  solute  Hamiltonian  Hs  according  to 


where  T is  the  kinetic  energy  operator  of  the  electrons  and  G is  the  Fock  potential. 
Modification  of  the  one  electron  part  of  the  solute  Hamiltonian  can  be  calculated 
as  the  penetration  of  the  point  charges  into  the  atomic  cores  of  the  solute  nuclei. 
It  is  this  method,  that  is,  the  change  of  Hs  due  to  the  penetration  integrals  within 
the  INDO  formulation  that  is  used  in  the  calculations,  the  results  of  which  are 
presented  in  the  next  section. 

4.3  Electronic  Spectra  of  a Solute  with  a Discrete  Number  of 
Solvent  Molecules 

The  solvent  induced  shifts  were  calculted  using  the  method  presented  in 
section  4.2.  The  origins  for  the  shifts  tabulated  below,  were  the  gas  phase 
values  calculated  from  INDO  and  those  obtained  experimentally.  Each  cluster 
was  simulated  as  described  in  the  introduction  to  this  chapter,  that  is,  by  building 
the  solvent  shell,  one  molecule  at  a time  till  the  interaction  energy  per  solvent 
molecule  converges  to  within  half  a kilocalorie  per  mole.  The  size  of  the  active 
space  in  the  configuration  interaction  procedure  was  restricted  to  ten  orbitals, 
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five  occupied  and  five  unoccupied,  the  same  space  used  in  all  the  calculations 
presented  in  Chapter  2. 

TABLE  4.1  Gas  phase  spectra  of  some  solute  molecules. 


SYSTEM 

CALCULATED 

EXPERIMENTAL 

TRANSITION  (cm'1) 

TRANSITION  (cm1) 

phenol 

36950 

36352  (ref  37) 

benzene 

37974 

38090  (ref  73) 

pyrazine 

29071 

30875  (ref  77) 

pyrimidine 

33235 

31073  (ref  24) 

The  first  system  investigated  was  the  step-wise  hydration  of  phenol  and  the 
results  are  given  in  Table  4.2. 

TABLE  4.2  Shifts  in  the  7r— >7r*  transitions  for  phenol  — (H2O),, 
n = 1,  2,  3,  5 


n 

CALCULATED  SHIFT 

EXPERIMENTAL 

(cm'1) 

SHIFT  (ref  37)(cm'1) 

1 

-123 

-356 

2 

-225 

-93 

3 

-196 

-89 

. 

• 

• 

5 

-99 

• 
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Experimental  results  using  the  MPS  techniques  (37)  show  a red  shift  in  the 
lowest  7T— +7 r*  band  of  phenol  when  hydrated  by  one  water  molecule.  This  is  due 
to  the  proton  donor  property  of  the  central  molecule  resulting  in  the  formation  of 
the  phenolate  ion.  The  first  water  molecule  is  therefore  expected  to  orient  itself 
such  that  the  phenolic  hydrogen  can  bond  to  the  oxygen  atom  of  the  water  mole- 
cule. Figure  4. 1 shows  the  structure  of  the  phenol-water  complex  as  calculated  by 
the  simulated  annealing  method.  The  arrangement  of  the  solvent  molecule  clearly 
shows  its  proton  acceptor  property. 

Further  hydration  by  another  molecule  causes  a blue  shift,  with  respect  to 
the  lowest  band  of  the  phenol-water  complex,  in  the  experimental  spectrum  for 
the  same  tt—ht*  band  as  listed  in  Table  4.2.  In  this  case  the  second  water  mol- 
ecule stabilises  the  phenolate  ion  resulting  in  the  observed  shift  as  verified  by 
Kasha  et  al  (72).  The  additional  solvent  molecule  should  orient  itself  such  that 
the  hydrogen  atoms  now  behave  as  proton  donors. 

The  calculated  spectral  shifts  show  that  two  water  molecules  are  necessary 
for  the  formation  of  the  phenolate  ion  which  is  then  stabilised  by  the  third  solvent 
molecule.  It  is  therefore  expected  that  two  water  molecules  will  behave  as  proton 
acceptors  towards  the  phenol.  Figures  4.2  and  4.3  show  the  phenol-(H20)2  and 
phenol-(H20)3  structures  calculated  by  the  annealing  method.  The  orientation  of 
the  first  two  water  molecules  in  figure  4.2  shows  their  acidic  behaviour.  In  figure 
4.3,  the  third  water  molecule  is  not  in  the  immediate  vicinity  of  the  phenolic 
group.  Spectroscopic  evidence  summarized  in  Table  4.2,  shows  that  the  spectral 
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shift  due  to  this  water  molecule  is  small,  and,  it  has  been  postulated  (37),  that  only 
two  water  molecules  are  closely  associated  with  the  phenol  molecule,  the  other 
water  begins  to  form  part  of  the  bulk  solvent,  a result  reproduced  by  the  annealing 
method.  The  blue  shift  in  the  spectrum  continues  with  increasing  hydration,as  the 
results  for  phenol-(H20)5  indicate. 

The  orientations  of  the  first  three  water  molecules  around  the  phenol  are 
given  in  Figures  4.1,  4.2  and  4.3. 


MPS  studies  of  the  interaction  of  benzene  and  diaza  compounds  with  small 
molecules  have  been  reported  recently  (74-77).  For  the  interaction  of  benzene 
with  one  molecule  of  water  and  ammonia,  the  calculated  interaction  energy  using 
the  annealing  method  and  shifts  in  the  7r  — ► 7r*spectrum  are  given  in  Table  4.3 
along  with  the  experimental  values  for  the  shifts  . The  geometry  of  the  ammonia 
molecule  was  taken  from  reference  (78). 


Figure  4.1  Orientation  of  H20  around  phenol. 
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Figure  4.2  Structure  of  phenol  - (H20)2- 


Figure  4.3  Arrangement  of  (H20)3  around  phenol. 
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TABLE  4.3  Calculated  and  experimental  shifts  for  benzene  with  H2O 
and  NH3. 


SYSTEM 

INTERACTION 

ENERGY 

(kcals/mole) 

CALCULATED 
SHIFT  (cm1) 

EXPERIMENTAL 
SHIFT  (cm1) 

Benzene  - H2O 

-2.66 

0 

78  (ref  75) 

Benzene  - NH3 

-4.49 

0 

-20  - +20  (ref  75) 

Figure  3.4  shows  the  structure  of  the  benzene-F^O  complex.  The  orientation 
of  the  ammonia  molecule  over  the  hydrocarbon  is  very  similar  to  the  former  case, 
a near  symmetrical  arrangement.  Thus,  no  change  in  the  electric  dipole  moment 
is  predicted  and  hence  no  spectral  shift  can  be  calculated. 

Pyrazine  interacting  with  NH3  and  [CH^Jn,  {n  = 1,  2,  3}  was  investigated 
next.  For  pyrazine/tCHUh  experimental  studies  (74)  show  two  n — > n*  peaks, 
one  corresponding  to  the  isotropic  structure  in  which  the  two  methane  molecules 
are  on  either  side  of  the  pyrazine  ring,  and  the  other  peak  due  to  the  anisotropic 
arrangement  of  the  methane  molecules.  The  simulated  annealing  method  could 
not  predict  the  isotropic  structure  for  this  system  although  it  is  expected  to  be  the 
more  stable  one.  The  shifts  in  the  n — ► n*  transition  peaks  and  the  calculated 
interaction  energies  are  tabulated  below. 
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TABLE  4.4  Interaction  energies  and  shifts  in  the  n— *7r*  transition  peaks  for 
pyrazine  with  NH3  and  methane. 


SYSTEM 

INTERACTION 

CALCULATED 

EXPERIMENTAL 

ENERGY 

SHIFT  (cm1) 

SHIFT  (cm1)  (ref 

(kcals/mole) 

76) 

pyrazine  - NH3 

-5.54 

-15 

117 

pyrazine  - CH4 

-1.19 

-7 

-32 

pyrazine  - (CH^ 

-2.86 

-9 

-38 

pyrazine  - (CH^ 

-5.01 

-21 

-70 

Figure  4.4  shows  the  position  of  the  NH3  molecule  with  respect  to  the  pyrazine 
ring. 


Figure  4.4  Position  of  NH3  above  the  pyrazine  ring. 
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Pyrimidine  interacting  with  NH3  and  methane  was  the  other  systems  studied. 
One  of  the  differences  between  this  system  and  the  pyrazine  / [CH^  is  that 
both,  the  isotropic  and  anisotropic  structures  for  its  interaction  with  methane  were 
calculated.  The  interaction  energies  and  the  shifts  in  the  n — > 7r*  transitions  are 
given  in  Table  4.5 

TABLE  4.5  Interaction  energies  and  shifts  in  the  n-*7r*  transition  peaks 
for  pyrimidine  with  NH3  and  methane. 


SYSTEM 

INTERACTION 

ENERGY 

(kcals/mole) 

CALCULATED 
SHIFT  (cm1) 

EXPERIMENTAL 
SHIFT  (cm1)  (ref 
74) 

pyrimidine  - NH3 

-6.13 

225 

285 

pyrimidine  - CH4 

-1.18 

33 

-57 

pyrimidine  - 
(CH4>2  (iso) 

-2.37 

75 

-112 

pyrimidine  - 
(CH4>2  (aniso) 

-1.75 

31 

-47 

The  anisotropic  and  isotropic  structures  of  pyrimidine  with  [CH^  are  shown 


in  Figures  4.5  and  4.6  respectively. 
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Figure  4.5  Anisotropic  structure  of  pyrimidine  - (CH^. 


Figure  4.6  Isotropic  structure  of  pyrimidine  - (CILjh. 

The  experimental  shifts  for  the  systems  studied  in  Tables  4.3,  4.4  and  4.5 
are  small,  and  can  be  due  to  geometry  changes  in  the  central  molecule  due  to  the 
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solvent  interactions.  Such  changes  have  been  ignored  in  the  calculations  presented 
here.  The  emphasis  in  calculating  these  shifts  is  to  investigate  the  power  of  the 
INDO  method  in  predicting  the  direction  of  the  shifts,  that  is,  towards  the  red  or 
blue  side  of  that  of  the  central  molecule. 

The  calculations  for  the  pyrazine-methane  clusters  reproduce  the  experimen- 
tal shifts  and  the  correct  trend  with  increasing  number  of  methane  molecules. 
However,  the  shift  calculated  for  the  pyrazine-ammonia  is  opposite  to  the  experi- 
mental observation.  In  the  case  of  pyrimidine  (Table  4.5),  the  trend  is  opposite  — 
the  shift  for  the  pyrimidine-NH3  complex  is  calculated  in  accordance  with  experi- 
ment whereas  the  pyrimidine-methane  complexes  are  predicted  to  have  blue  shifts. 
Increasing  the  size  of  the  configuration  interaction  space  to  twenty  molecular  or- 
bitals, ten  occupied  and  ten  unoccupied,  produced  the  same  trend  in  the  results. 

The  method  fails  for  clusters  in  which  the  solute/solvent  systems  comprise 
polar/nonpolar  molcular  combinations.  The  hydrogen  bonding  capabilities  of  di- 
azines  in  their  excited  states  has  been  experimentally  verified  (75,79).  In  cases 
like  pyrazine  with  NH3,  the  method  used  to  calculate  the  spectral  shifts  should 
be  able  to  describe  such  intermolecular  hydrogen  bonding,  which  stabilises  the 
excited  state.  The  extent  of  the  stabilisation  may  be  small,  but  is  indeed  important 
in  the  calculations  as  the  shifts  measured  are  less  than  250  wavenumbers.  Re- 
placing the  solvent  molecule  by  point  charges  removes  any  molecular  properties 
the  solvent  may  have,  which  in  turn  neglects  specific  intermolecular  effects  (76). 
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4.4  Shifts  in  the  Electronic  Spectra  of  a 


Solute  Using  the  Discrete-Continuum  Model 


Using  equation  (4.1)  and  the  continuum  method  presented  in  Chapter  2,  the 
Hamiltonian  H’  for  the  solute  - discrete  solvent  complex,  centered  in  a cavity  of 
radius  ‘ b ’ in  the  continuum  solvent  is  now  given  by 


In  order  to  treat  the  solute  - discrete  solvent  cluster  within  the  continuum 
model,  a value  for  the  radius  of  the  cavity  which  contains  the  cluster  is  required. 
Two  methods  were  used  to  evaluate  this  quantity  — one  was  to  calculate  the 
radius  rMD  as  1.5  times  the  radius  obtained  from  the  mass  density  as  explained  in 
Chapter  2.  The  other  method  was  to  take  the  maximum  and  minimum  distances, 
rmax  and  rmin,  of  the  center  of  mass  of  the  solvent  molecules  from  the  origin. 
Once  the  value  for  the  radius  has  been  determined,  the  spectra  can  be  calculated 
using  the  Hamiltonian  given  in  equation  (4.4)  and  the  continuum  solvent  model 
detailed  in  Chapter  2. 

The  calculated  and  experimental  values  for  the  shifts  in  the  spectral  transitions 
for  pyridine  /[H20]n  {n  = 1,  2,  4 and  6}  for  the  various  values  for  the  cavity 
radius  are  given  in  Table  4.6. 


N 


N 


4.4 
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TABLE  4.6  Cavity  radii  and  shifts  in  the  n—>n*  spectra  for  pyridine  / (H20)n 
n = 1,  2,  4 and  6 


SYSTEM 

CAVITY  RADIUS 
(A) 

CALCULATED 
SHIFT  (cm1) 

EXPERIMENTAL 
SHIFT  (cm1) 

pyridine/H20 

tmd  = 4.76 

-94 

I’m  ax  = Imin~  3.78 

2830 

pyridine/(H20)2 

fMD  = 4.76 

-147 

rmax=  4.39 

512 

Tmin  = 4.20 

219 

pyridine  /(H2CO4 

tmd  = 4.76 

-125 

fmax=  6.26 

139 

rmin  = 5.38 

-317 

pyridine/(H20)6 

fMD  = 4.76 

346 

fmax  = 16.06 

4 

rmin=  12.76 

9 

pyridine/water 

713  (ref  28) 

The  results  of  a similar  treatment  for  phenol/(H20)n  {n  = 1,  2,  3}  are  given 
in  Table  4.7.  The  calculated  shifts  in  the  spectral  transitions  for  the  phenol/H20 
system  using  the  distance  from  the  solute  origin  to  the  center  of  mass  of  the 
solvent  molecule  as  the  cavity  radius  is  ridiculous,  the  i r*  state  is  predicted  to 
have  a dipole  moment  of  224  D.  A similar  result  is  obtained  for  the  phenol/ 
(H20)2  system.  It  was  therefore  decided  to  use  the  scaled  radius  obtained  from 
mass  density  measurements  as  explained  in  Chapter  2. 
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TABLE  4.7  Shifts  in  the  n — > 7r*  spectrum  for  phenol  / (ThO),, 
n = 1,  2,  3,  5 


SYSTEM 

CALCULATED  SHIFT 
(cm’1) 

EXPERIMENTAL  SHIFT 
(cm'1) 

phenol/H20 

-66 

phenol  / (H20>2 

-192 

phenol/(H20)3 

-152 

phenol/(H20)5 

29 

phenol/water 

685  (ref  33) 

From  the  results  in  Tables  4.6,  the  best  value  for  the  cavity  radius  is  that 
obtained  from  mass  density  measurements.  It  was  therefore  decided  to  use  this 
values  for  each  of  the  solutes  studied  using  this  model. 

In  the  case  of  pyridine  with  hexane  and  ethanol,  the  cluster  buildup  was 
performed  using  the  simulated  annealing  method.  In  each  instance,  one  solvent 
molecule  was  found  to  satisfy  the  criterion  for  the  convergence  of  the  interaction 
energy  of  the  cluster. 

TABLE  4.8  Shifts  in  the  tt  — ► n*  spectral  transitions  using  the 
discrete-continuum  approach. 


SYSTEM 

CALCULATED 

SHIFT(cm'!) 

EXPERIMENTAL  SHIFT 
(cm'1) 

pyridine/water 

346 

713  (ref  28) 

pyridine/hexane 

24 

1491  (ref  30) 

pyridine/ethanol 

307 

560  (ref  29) 

phenol/water 

29 

685  (ref  33) 
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4.5  Conclusions 


Simulation  of  solvent  clusters  around  the  nitrophenols  could  not  be  carried 
out  as  no  parameters  are  available  for  the  nitrogen  atom  in  the  nitro  group. 

The  agreement  between  the  calculated  and  experimental  shifts  presented  in 
Table  4.8  is  poor,  although  the  discrete  - continuum  model  is  expected  to  be  a 
better  description  of  the  solvent  than  the  continuum  model  alone.  Nevertheless, 
the  calculated  shift  for  pyridine  in  water  has  improved  compared  to  the  value 
of  290  wavenumbers  obtained  for  the  continuum  model.  For  the  other  cases, 
however,  there  is  no  improvement. 

In  these  calculations,  the  solute  molecule  was  considered  to  be  a monomer, 
an  assumption  that  may  not  be  valid  for  solutes  where  strong  associations  such 
as  hydrogen  bonding  are  possible.  Phenol  is  an  example  of  a molecule  that 
exists  as  dimers  and  trimers  in  water  (37).  Furthermore,  phenol  which  exhibits 
acidic  behaviour  in  water,  exists  as  a mixture  of  phenol  and  the  phenolate  ion, 
an  effect  not  considered  in  the  calculations.  Experimental  values  for  these  shifts 
are  a source  of  controversy  in  the  literature  as  the  purity  of  the  solvent  water  is 
extremely  important. 

As  explained  in  Chapter  2,  the  calculated  spectral  shifts  are  dependent  on  the 
value  of  the  cavity  radius.  The  method  that  was  used  to  evaluate  this  parameter 
is  the  same  for  all  the  solute  molecules.  Such  an  approach  may  not  be  useful  and 
it  is  therefore  necessary  to  be  able  to  define  a cavity  for  the  solute  which  has  a 
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more  realistic  shape  and  suitable  dimensions.  In  the  following  chapter,  such  an 
algorithm  and  the  results  of  its  applications  are  presented. 


CHAPTER  5 

AN  ALGORITHM  FOR  A MORE  REALISTIC  DESCRIPTION 
OF  THE  SOLUTE  MOLECULE 


5.1  Introduction 

In  Chapters  2 and  4,  the  solute  molecule  was  assumed  to  be  a point  dipole 
centered  in  a spherical  cavity  having  a suitable  radius.  Very  few  molecules  can  be 
considered  to  be  even  approximately  spherical  in  shape  and,  if  its  structure  has  side 
chains,  the  spherical  cavity  model  is  certainly  a poor  description  of  the  solute. 
Furthermore,  in  studying  processes  like  molecular  dissociation,  it  is  unrealistic 
to  have  the  molecular  fragments  in  a single  spherical  cavity  surrounded  by  the 
solvent. 

It  is  therefore  necessary  to  determine  a more  realistic  shape  for  the  molecular 
cavity,  and  the  algorithm  that  has  been  used  in  this  chapter  is  based  on  the 
method  presented  by  Miertus  et  al  (80,  81).  The  molecular  systems  studied 
are  the  dissociation  of  HF  and  the  hydration  of  the  neurotransmitter  7 - amino 
butyric  acid  (GABA). 

5.2  Mathematical  Formulation  of  the  Algorithm 

According  to  the  method  of  Miertus  et  al,  every  atom  in  the  solute  molecule 
is  centered  in  a sphere  having  the  van  der  Waals  radius  of  that  particular  atom. 
The  cavity  that  surrounds  the  solute  molecule  is  thus  approximated  as  a set  of 
interlocking  spheres. 
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Suppose  the  solute  molecule  has  a gas  phase  charge  distribution  p°  corre- 
sponding to  the  Hamiltonian  H0.  When  this  solute  is  immersed  in  a solvent 
having  a dielectric  constant  e,  po  would  induce  a set  of  image  charges  { qi°  } , i = 
1,  2,  N,  on  the  surfaces  of  the  interlocking  spheres  which  define  the  boundary 
of  the  cavity.  The  value  of  these  charges  will  depend  on  the  value  of  the  dielectric 
constant  — higher  the  value  of  e,  larger  the  magnitude  of  the  charge.  The  number 
of  such  charges,  N,  is  predetermined  by  the  user. 

These  image  charges  will  polarise  each  other  resulting  in  a new  set  of  charges 
{ qj'  } and  the  solute  molecule  will  experience  the  effect  of  these  modified  charges. 
Thus,  the  new  Hamiltonian  for  the  solute  is  H where 

H'  = H0  + V 5.1 

where  V is  the  effect  of  the  image  charges.  The  new  charge  distribution  p 
corresponding  to  H'  will  now  induce  image  charges  { q,"  } which  effects  H'  in  a 
manner  very  similar  to  that  of  p°  on  Ho-  This  process  of  modifying  the  gas  phase 
Hamiltonian  of  the  solute  using  the  influence  of  the  induced  image  charges,  V, 
is  repeated  until  self-consistency  in  the  total  energy  of  the  solute  with  respect  to 
the  Hamiltonian  H*  is  obtained. 

As  a first  step  in  applying  this  algorithm,  an  initial  set  of  point  charges 
{ qj°  } i = 1,  2,  ....  N , is  required.  These  charges  which  are  located  on  the 
surfaces  of  the  spheres  are  associated  with  the  net  charge  qA  and  the  van  der 
Waals  radius  rA  of  the  particular  atom  A on  which  the  sphere  is  centered.  The 
surface  charge  density  <7j°  at  the  center  of  the  surface  element  j on  the  sphere 
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centered  on  atom  A is  given  by  (82) 

= ~{J^u/Na  5.2 

erA 

where  NA  is  the  number  of  point  charges  around  atom  A.  The  charges  qj°  are 
now  given  by 

q°  = a)  A Sj  5.3 

where  ASj  is  the  element  of  surface  area  associated  with  qj°  given  by 
AS]  = 2 r\  sin  9aj  sin(A0/2)  A <f>  5.4 


#A,j  is  the  angle  made  by  the  radius  of  the  sphere  to  the  center  of  the  surface 
element  under  consideration  with  the  positive  z-axis,  and  A#  and  A<f>  are  the 
increments  in  the  spherical  angular  coordinates  that  describe  the  sphere.  In  order 
that  the  surface  elements  have  the  same  area  at  axial  and  equatorial  positions, 
the  increment  in  the  angle  <f>,  A <f>,  is  calculated  as  a function  of  the  angle  9 and 
A 9 according  to 


Ay>  = 


A 9 
sin  9 


Substituting  equation  (5.5)  in  (5.4)  yields 


5.5 


AS]  = 2 r\  sin  ( A0/2)  A 9 5.6 


Thus,  for  a given  increment  of  A0,  equal  areas  are  evaluated  on  the  surface  of 
the  van  der  Waals  spheres  and  at  the  center  of  each  of  these  elements  is  located 
a point  charge  qj°  calculated  according  to  equations  (5.2)  and  (5.3). 
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The  polarisation  of  each  image  charge  qj°  in  the  presence  of  all  the  charges 


is  evaluted  next  as  the  sum  of  two  effects  — one  due  to  the  presence  of  all 
the  other  point  charges  except  itself  and  the  other  contribution  attributed  to  the 
self-polarisation  of  qj°  whose  effect  is  to  reduce  the  magnitude  of  the  charge. 

The  first  effect  is  calculated  using  a finite  difference  method.  The  gradient  of 
the  potential  between  point  charges  qj°  and  qj°,  i not  being  equal  to  j,  is  evaluated 
as  the  difference  in  the  potential  at  points  r,  and  rj  + 6,  where 
I 6 I was  arbitrarily  chosen  to  be  0.01  A.  The  self-polarisation  term  was  then 
evaluated  taking  into  account  the  curvature  of  the  sphere  and,  the  process  of 
charge  polarisation  is  repeated  till  the  charge  distribution  has  converged  within  a 
predetermined  threshold  value.  The  equation  used  is 


rj  and  Fi  are  the  position  vectors  of  the  image  charges  j and  i respectively  and 
m is  the  number  of  cycles  required  for  convergence. 

The  effect  of  the  point  charges,  V,  which  modifies  the  Hamiltonian  Ho 
according  to  equation  (5.1)  is  defined  by 


The  summation  index  k runs  over  all  atomic  cores  of  charge  Zk  and  j represents 
the  electrons  in  the  solute,  {q,*}  are  the  polarised  image  charges  on  the  surface 
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of  the  Van  der  Waals  spheres.  The  effect  of  V,  as  defined  in  the  above  equation, 
is  to  modify  the  one  electron  part  of  the  solute  Hamiltonian  Ho- 

In  defining  the  Van  der  Waals  spheres  around  each  atom,  the  tail  of  the  charge 
distribution  p( r)  of  the  solute  molecule  will  lie  outside  the  boundary  of  the  cavity. 
As  a result  of  this  truncation,  the  sum  of  all  the  image  charges  { q**  } will  not  be 
zero.  It  is  therefore  necessary  to  normalise  these  point  charges.  We  choose 
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where  Q is  the  sum  of  the  absolute  values  of  the  image  charges  and  qi,+  and  qi,. 
represent  the  positive  and  negative  charges  respectively. 

The  method  outlined  above  is  incorporated  into  the  INDO  method  (12,  13) 
and  the  results  obtained  are  presented  in  the  next  section. 


5.3  Results  and  Discussion  of  the  Applications  of  the  Point 

Charge  Model 

The  first  system  studied  was  the  dissociation  of  HF  in  water.  Gas  phase 
dissociation  favours  the  homolytic  cleavage  of  the  HF  bond  (83)  resulting  in 
the  formation  of  H and  F radicals  while  heterolytic  fission  is  observed  in  polar 
solvents.  Previous  work  on  this  system  (20)  used  the  reaction  field  theory  of 


93 


Onsager  in  which  both  atoms  are  in  the  same  spherical  cavity  as  the  bond  distance 
increased.  The  radius  of  the  cavity  did  not  change  as  a function  of  the  distance 
between  the  H and  F atoms. 

In  the  calculations  presented  here,  the  total  number  of  point  charges  on  the 

cavity,  N,  and  the  increment  A 0 in  the  spherical  angle  6 is  varied.  The  energy 

E ( in  a.u.  ) and  the  formal  charges  on  the  atoms  are  evaluated  as  a function  of 

the  interatomic  distance.  In  Table  5.1,  the  energy  as  a function  of  the  interatomic 

distance  is  tabulated.  Table  5.2  gives  the  formal  charge  on  the  hydrogen  atom, 

6h,  as  a function  of  the  bond  distance  and  Figure  5.1  shows  the  potential  energy 

curve  for  the  dissociation  of  HF  in  water. 

TABLE  5.1  Energy  of  HF  (a.u.)  as  a function  of  the  interatomic  distance  r,  N 
and  increment  in  the  spherical  angle. 


THF  (A) 

GAS 

N = 10, 

N = 100, 

N = 200, 

N = 240, 

PHASE 

o 

rH 

II 

<1 

A 9 = 10 

> 

II 

L/l 

II 

<1 

0.5 

-24.23971 

-24.29171 

-24.40536 

-24.47807 

-24.45333 

0.7 

-25.51578 

-26.32725 

0.8 

-25.74363 

-25.79341 

-26.22123 

-26.42758 

-26.32314 

1.0 

-25.85487 

-25.90582 

-26.17718 

-26.10647 

-26.08222 

1.2 

25.77653 

-26.00378 

1.5 

-25.59443 

-25.71323 

-25.99364 

-25.79994 

-25.82210 

2.0 

-25.37826 

-25.59435 

-25.77315 

-25.71218 

-25.73283 

2.5 

-25.27624 

-25.52577 

-25.58584 

-25.63215 

-25.65395 

3.0 

-25.22995 

-25.72313 

-25.84482 

-25.783427 

-25.80822 

4.0 

-25.19369 

-25.89520 

-25.93880 

-25.94534 

-25.96936 
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TABLE  5.2  Formal  charge  on  the  hydrogen  atom  as  a function  of  the 
bond  distance  r. 


rHF  (A) 

N = 10,  AO  = 

N = 100,  AO 

N = 200,  A0 

N = 240,  A0 

10 

= 10 

= 5 

= 5 

0.5 

0.744 

0.808 

0.803 

0.799 

0.8 

0.479 

0.508 

0.573 

0.564 

1.0 

0.421 

0.522 

0.504 

0.499 

1.5 

0.541 

0.753 

0.604 

0.625 

2.0 

0.733 

0.884 

0.843 

0.861 

2.5 

0.822 

0.891 

0.929 

0.943 

3.0 

0.996 

0.999 

0.998 

0.998 

4.0 

1.000 

1.000 

1.000 

1.000 

The  results  in  Table  5.2  clearly  show  that  the  strong  dielectric  nature  of  the 
solvent  favors  the  formation  of  the  proton  and  fluoride  ion.  The  accumulation 
of  charge  on  the  atoms  increases  dramatically  for  values  of  the  bond  distance 
greater  than  1.5  A.  When  the  interatomic  distance  is  2.5  A,  the  molecule  is  almost 
completely  dissociated,  and  at  4A,  the  proton  and  the  flouride  ion  are  the  species 
found  in  aqueous  solution.  The  increase  in  the  energy  when  the  bond  distance 
is  about  2.5  A is  due  to  the  heat  of  hydration,  an  effect  observed  by  Karelson 
et  al  (20). 
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The  neurotransmitter  7-aminobutyric  acid  (GABA)  was  studied  next  In  polar 
solvents  the  zwitterionic  structure  H3N+(CH2)3COO~  is  observed  as  the  solvent 
can  stabilise  the  charge  separated  structure.  The  charge  densities  on  the  carboxylic 
carbon  atom,  denoted  />(Ci)  and  nitrogen  atom,  p(N)  were  calculated  for  the 
zwitterion  as  a function  of  the  number  of  image  points  N,  and  increment  A 8 in 
the  spherical  angle  9 using  the  method  presented  in  this  chapter.  The  results  are 
given  in  Table  5.3 


96 


TABLE  5.3  Charge  densities  on  the  carboxylic  carbon  and 
nitrogen  atoms  of  GABA. 


GAS 

N = 16, 

N = 32, 

N = 48, 

N = 64, 

N = 160, 

PHASE 

A e = 30 

AS  = 10 

Ad  = 10 

A 6 = 10 

A 6 = 10 

Energy 

-75.9005 

-76.5136 

-77.2897 

-77.1395 

-77.1946 

-77.1384 

(a.u.) 

p(Ci) 

0.444 

0.408 

0.352 

0.363 

0.380 

0.391 

P( N) 

-0.317 

-0.282 

-0.289 

-0.292 

-0.294 

-0.293 

The  positive  nature  of  the  carboxylic  carbon  is  reduced  as  the  number  of  image 
charge  points  increase  as  is  the  negative  charge  density  on  the  nitrogen  atom. 

The  energies  of  GABA  and  its  zwitterion  in  the  gas  phase  and  water  calculated 
using  the  algorithm  presented  in  this  chapter  are  given  below. 

TABLE  5.4  Energies  (a.u.)  of  GABA  and  its  zwitterion. 


GAS  PHASE 

SOLVATED  N=160,  A 6 
= 5 

GABA 

-76.0419 

-76.2645 

zwitterionic  form 

-75.9004 

-77.1383 
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The  above  results  indicate  that  the  zwitterionic  structure  is  more  stable  in 
water  than  GABA  itself  in  accordance  with  experiment.  This  is  to  be  expected 
as  water  being  a polar  solvent  can  stabilise  the  charge  separated  structure.  In  gas 
phase  ths  observation  is  reversed  as  no  medium  is  available  for  the  stabilisation 
of  the  charges  in  the  zwitterion. 

This  result  is  certainly  gratifying  as  a method  as  simple  as  that  presented  in 
this  chapter  is  able  to  reproduce  experimental  observations.  In  a previous  study  if 
GABA  and  its  zwitterion  using  the  Onsager  cavity  model  (84),  the  reaction  field 
factor  had  to  be  manipulated  in  order  to  get  results  which  agree  with  experiment. 

The  method  presented  in  this  chapter  will  not  be  able  to  show  the  formation 
of  the  zwitterionic  structure  as  the  explicit  nature  of  the  solvent  molecules  and 
the  electronic  processes  involved  in  intermolecular  interactions  will  have  to  be 
considered  explicitly.  Such  a treatment  requires  supermolecular  calculations  and 
have  been  reported  by  Purvis  and  Zemer(85). 

An  algorithm  which  divides  the  surface  area  of  the  sphere  into  equal  segments 
and  determines  the  increase  in  the  spherical  angle  <j>  as  a function  of  8 is  one  way 
in  which  the  method  presented  here  can  be  improved. 


CHAPTER  6 
CONCLUSIONS 


The  results  of  the  calculations  presented  in  this  thesis  have  been  discussed  in 
detail  in  the  relevant  sections  in  the  proceeding  chapters.  Only  general  comments 
on  the  applicability  of  the  methods  and  the  theory  used  will  be  presented  here. 

In  spite  of  its  inherent  assumptions,  the  reaction  field  method  has  been  shown 
to  be  useful  in  obtaining  solvent  induced  spectral  shifts  for  a series  of  molecules. 
In  those  cases  where  strong  interactions  between  the  solute  and  the  solvent  ex- 
ist, for  example,  pyridine  in  water,  the  method  is  not  expected  to  be  successful. 
The  theory  on  which  the  method  is  based  , that  of  weak  solute-solvent  interac- 
tions, neglects  any  change  in  the  chemical  structure  of  the  constituent  systems. 
Thus,  phenomena,  like  the  formation  of  a pyridinium  ion,  cannot  be  adequately 
described  using  the  reaction  field  method. 

A drawback  in  the  method  is  that  the  results  are  dependent  on  the  value 
of  the  cavity  radius,  as  the  solvent  perturbation,  that  is,  the  reaction  field 
term  is  inversely  proportional  to  the  cube  of  this  parameter.  If  the  radius 
is  too  small,  the  effect  of  the  solvent  is  overestimated;  too  large  a value 
will  not  include  the  solvent  effects  in  any  meaningful  manner,  thereby  re- 
producing the  vacuum  state  properties  of  the  solute.  It  is  reasonable  to  ex- 
pect that  any  algorithm  that  is  used  to  calculate  the  value  for  the  cavity 
radius  should  include  some  phsical  property  of  the  solute.  However,  re- 
cent studies  using  a multiconfigurational  selfconsistent  reaction  field  technique 
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(86)  have  shown  that  properties,  like  the  ionisation  potential  of  a solvated 
molecule,  can  be  evaluated  with  resonable  accuracy  using  the  reaction  field 
scheme. 

The  continuum  solvent  model  assumes  the  solute  to  be  a monomer  at  infinite 
dilution.  Experimentally  this  may  not  be  the  case  as  the  entity  whose  spectra 
is  being  measured  could  be  an  associated  form  of  the  solute,  an  effect  which  is 
important  in  dealing  with  solutes  with  hydrogen  bonding  capabilities. 

The  success  of  the  simulated  annealing  method  using  the  Fraga  intermolecular 
potential  energy  function  and  the  limitations  in  its  applicability  have  been  demon- 
strated in  Chapter  3.  As  the  annealing  parameter  T*  has  units  of  energy  and  not 
temperature,  the  simulated  annealing  method  is  well  suited  for  incorporation  into 
standard  quantum  chemical  methods  used  in  evaluating  molecular  properties. 

From  a computational  point  of  view,  the  time  required  to  evaluate  minima  is 
somewhat  longer  than  in  a gradient  search  technique,  as  the  annealing  algorithm 
has  to  be  given  an  adequate  number  of  moves  and  cycles  in  order  to  ensure  that 
the  potential  energy  surface  has  been  scanned  as  completely  as  possible.  How- 
ever, the  advantage  is  that,  given  these  conditions,  the  global  minimum  is  reached 
as  explained  in  Chapter  3. 

The  discrete-continuum  model  for  the  solvent  includes  the  anisotropy  of  the 
solvation  shells  and  the  bulk  effects  of  the  solvent.  The  formation  of  the  solvent 
cluster  is  time  consuming,  especially  for  larger  solvent  molecules  with  conforma- 
tional flexibilty  as  in  the  case  of  ethanol  and  hexane.  In  the  calculations  presented 
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in  chapter  4 for  pyridine  in  ethanol  and  hexane  such  conformational  effects  were 
not  considered. 

Recent  spectroscopic  studies  have  been  reported  for  the  stepwise  solvation  of 
tryptophan  and  other  indole  derivatives  (87).  The  formation  of  the  solvent  shells 
would  favour  a buildup  around  the  saturated  side  chains  and  therefore,  the  effect 
on  the  7r  sytem  will  be  very  small.  The  method,  as  presented  in  this  dissertation 
is  not  expected  to  be  able  to  evalute  the  spectral  shifts  in  such  cases  with  any 
reasonable  accuracy. 

Some  improvements  to  the  methods  will  now  be  presented.  A better  algorithm 
for  the  calculation  of  the  cavity  radius  is  required.  The  method  should  not  be 
time  consuming,  and,  furthermore,  need  not  be  solely  dependent  on  the  solute 
properties.  Inclusion  of  some  of  some  solvent  parameters  like  the  surface  tension 
seems  reasonable  from  physical  considerations.  The  dispersion  term  and  the 
effects  of  the  polarisation  of  the  solute  by  the  solvent  will  have  to  be  considered 
explicitly  so  as  to  make  the  reaction  field  method  applicable  to  both  polar  and 
nonpolar  systems. 

The  success  of  the  intermolecular  potential  function  used  here  in  treating 
neutral  systems  has  been  previously  emphasised.  However,  in  order  to  study 
charged  systems  and  molecules  with  a large  separation  of  charges,  like  zwitterions, 
another  function  which  explicitly  considers  the  interactions  in  such  systems  will 
have  to  be  used. 

The  defining  of  a more  realistic  shape  for  the  solute  cavity  in  terms  of 
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interlocking  van  der  Waals  spheres  makes  it  possible  to  study  a wider  range  of 
molecules,  specially  those  with  side  chains.  The  results  obtained  will  be  dependent 
on  the  values  of  the  radii  used  for  these  spheres. 

The  results  presented  in  this  thesis  demonstrate  that  it  is  indeed  possible  to 
construct  simple  and  useful  models  for  solvents  without  resorting  to  expensive 
supermolecular  calculations  whose  applicability  is  severely  restricted  by  the  size 
of  the  constituent  molecules.  The  incorporation  of  the  solvent  effects  into  the 
Hamiltonian  for  the  solute  can  be  done  in  various  physically  meaningful  ways 
and  this  flexibility  is  very  helpful  in  studying  problems  dealing  with  solvation. 


APPENDIX 

THE  REACTION  FIELD  METHOD 


In  this  method,  the  solvent  molcules  are  replaced  by  a continuum  and  the 
solute  is  treated  as  a polarisable  dipole  centered  in  a cavity  which  is  embedded  in 
the  continuum.  For  purposes  of  simplicity,  the  cavity  is  assumed  to  be  a sphere 
of  radius ' a As  the  solute  is  neutral,  the  potential  in  the  cavity  can  be  obtained 
by  solving  Laplace’s  equation.  Without  loss  of  generality,  the  z-axis  is  taken  to 
be  the  direction  of  the  dipole,  the  origin  of  the  system  of  axes  being  coincident 
with  the  center  O of  the  sphere. 

As  the  system  has  spherical  symmetry,  the  <f>  dependence  of  the  potential 
will  not  appear  explicitly  in  the  expression  for  the  potential  V.  Thus  V can  be 
regarded  as  a function  of  the  variables  r and  6. 

Let  Vj  and  V2  be  the  potentials  outside  and  inside  the  cavity  respectively.  As 
the  solute  and  solvent  molecules  are  neutral,  there  are  no  excess  charges  present. 
Therefore  Vj  and  V2  satisfy  the  Laplace  equation 


Suppose  a general  potential  V(  r ,6  ) which  satisfies  the  Laplace  equation  is 
separable  according  to 


A Vi  = 0 


A.  1 


AV2  = 0 


A.  2 


V{r,9)  = R(r)Q(6) 


A.Z 


The  Laplace  equation  is  given  in  spherical  polar  coordinates  by 


AA 
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Substituting  equation  (A.3)  in  (A.4)  followed  by  separation  of  variables  yields 


r d { dR{r) 


dr 


R (r)  dr 

where  A is  the  constant  of  separation. 

The  0 part  of  equation  (A.5)  is 

d2S 
de 2 


i d2e 

’0(0)  dd2 


= A A.5 


= — A0 


whose  solution  is  given  by 

0 = A cos  (%/jw)  + B sin  (VX»)  A.  6 

where  A and  B are  constants.  Using  the  periodic  properties  of  the  function  0(0) 
in  equation  (A.6)  yields 


A = n2  n = 0,  1,  2,  . . . 


Thus, 


Q — A cos  (n8)  + B sin  (n6)  A.l 


The  function  0(0)  can  be  written  in  terms  of  the  Legendre  polynomials  in  Cos(0), 
P n(Cos(9)). 

The  radial  part  of  equation  (A.5)  now  reads 


- n2R  (r)  = 0 

drz  dr 

which  is  an  Euler  equation  and  its  solution  is  given  by 


A.8 


R (r)  = Cir"  + C2r  " 


A.  9 
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Inserting  equations  (A.8)  and  (A.9)  in  (A.3)  and  using  the  principle  of  linear 
superposition,  one  obtains 


V = 


Pn  (cos  8 ) 


A. 10 


The  constants  Ln  and  Mn  incorporate  A,  B,  Ci  and  C2- 


Solution  for  the  potential  outside  and  inside  the  cavity. 

The  potentials  Vj  and  V2  outside  and  inside  the  cavity  are  given  by: 

°°  / A/f^X 

V,(r)=£  UL‘)(r)  + -^r  P.  (cos 0)  All 

n=0  V / 

°o  / ,,(2)\ 

Vj(r)  = £ A2)(r)  + -£y  P.(«»0)  A12 

n=0  V / 

In  order  to  determine  the  coefficients  Ln®,  Mn®  , { i = 1,  2 },  boundary  conditions 
on  the  potential  are  required.  These  conditions  are: 

(r)]^  = 0 A.13 

[Vi(r)]r=0  = [V2(r)]r=0  A.14 


'dVi(r)' 

'dV2(rf 

dr 

r=a 

dr 

where  e is  the  dielectric  constant  of  the  bulk  fluid.  The  only  source  of  charge 
inside  the  cavity  is  the  dipole  p.  Hence 


V2  (r)  = cos  8 
rz 


A.  16 
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where  6 is  the  angle  an  arbitrary  point  in  the  cavity  makes  with  the  positive 
z-axis.  Thus, 


mJ2)  = n 


AA1 


and 


Mi2)  = 0 Vn,  (n  ± 2)  A.  18 
Inserting  equations  A.17  and  A.18  into  A. 12  yields 

oo 

v2  (r)  = ^2  L^rnpn  (cos  0)  + -^  cos  6 ^AS 

n= 0 

From  equation  A.  13,  it  is  clear  that  in  order  to  have  a finite  potential, 

i!,1*  = 0 Vn 


This  simplifies  Vi(r)  as 

v’>(r)  = ESrp”<cos#)  A2° 

n= 0 

The  Legendre  polynomials  are  linearly  independent.  Using  equations  A.  14  and 
A.15  in  A. 19  and  A.20  yields 


n/1 


Mi1}  _ (2) 

a"+1  Ln 
and 

(2) 

- e (n  + 1)  = nl;  a 

Equations  A.21  are  consistent  if  and  only  if 

M{n]  = 0 n/1 


A.21 


L{n2)  = 0 
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For  the  special  case  when  n = 1,  equations  A.  19  and  A.20  together  with  the 
boundary  conditions  yield 


M. 


(i) 


= L^a  + ^ 

or 


A.  22 


and 


A. 

or  or 


23 


Li  (2)  and  Mi  (1)  can  now  be  obtained  from  these  two  equations  as 


and 


(2)  _ 2 (e  — 1)  /x 

1 (2e  + l)a3 


= STTT" 


A. 24 


A. 25 


Thus,  the  potentials  outside  and  inside  the  cavity  are  given  by 

Vl(r)  = —?—^cos6  A. 26 
v ' 2e  + 1 r2 


and 

V2  (r)  = 4 cos  0 - 2i£~!--4 r cos  9 A.27 

’ r 2 2e  + l a3 

From  equation  (A.27)  , it  is  clear  that  the  field  inside  the  cavity  is  enhanced  by 
a uniform  field  R given  by 


R = 


1 2(e  — 1) 
a3  2e  + l ^ 


A. 28 


This  field  is  called  the  reaction  field  and  is  due  to  the  the  image  created  by  the 
solute  dipole  in  the  surrounding  continuum  fluid. 
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